{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "33a0ad43-755a-47cf-b2af-fdd5fbe451e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from spn.structure.leaves.parametric.Parametric import Categorical\n",
    "\n",
    "spn = 0.4 * (Categorical(p=[0.2, 0.8], scope=0) *\n",
    "             (0.3 * (Categorical(p=[0.3, 0.7], scope=1) *\n",
    "                     Categorical(p=[0.4, 0.6], scope=2))\n",
    "            + 0.7 * (Categorical(p=[0.5, 0.5], scope=1) *\n",
    "                     Categorical(p=[0.6, 0.4], scope=2)))) \\\n",
    "    + 0.6 * (Categorical(p=[0.2, 0.8], scope=0) *\n",
    "             Categorical(p=[0.3, 0.7], scope=1) *\n",
    "             Categorical(p=[0.4, 0.6], scope=2))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b39e4565-e31d-46e1-952a-f12fa47c3ce3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from spn.algorithms.Marginalization import marginalize\n",
    "\n",
    "spn_marg = marginalize(spn, [1,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "144d646d-61ee-4501-8e68-fa5ed3749e5e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/80lEQVR4nO3de1SU95kH8O8MN7nNoAIionLVGcYr11GbCIpJWpvUpAkk0YTYk5ptQzY99XJ2u3vSbfbs9lTatFa33aNNxahpQtrENGmbRhRyEbkJGhgGEGZABLkpMFxnGObdP9yhXlDeGd6Z9zLP55w9pwvO+z48fTrPe/ldZAzDMCCEEEIERs53AIQQQsh0qEERQggRJGpQhBBCBIkaFCGEEEGiBkUIIUSQqEERQggRJGpQhBBCBIkaFCGEEEGiBkUIIUSQqEERQggRJGpQhBBCBIkaFCGEEEGiBkUIIUSQqEERQggRJGpQhBBCBIkaFCGEEEHy5jsAIkxWqxUGgwF6vR4NDQ0YGBjAxMQEfHx8EBISApVKBbVajdjYWHh7UxkRQrgnox11iR3DMKirq0NhYSGKi4sRGhqK2NhYREdHQ6lUwtvbG1arFYODg2htbYXBYEBfXx82bdqE7OxsaDQayGQyvv8MQohEUIMiAIDS0lIcPHgQw8PDyMrKQkZGBoKDg2f83NDQEIqLi1FUVASFQoG8vDysX7/eDRETQqSOGpSHGxoawhtvvIGKigrk5uYiKSkJcrnjryZtNhuqq6tRUFCA9PR07N69G0FBQS6ImBDiKWiQhAfT6XTIycmB2WzG/v37kZKS4lRzAgC5XI6UlBTk5+fDbDYjOzsbOp2O44gJIZ6E7qA81IULF7Bv3z7s2rULqampnB+/srIShw8fxv79+5GcnMz58Qkh0kcNygPpdDq8+uqreOWVV7By5UqXnae2thYHDx7EgQMHoNFoXHYeQog0UYPyMMPDw8jOzkZubq5L7pzuVFlZiWPHjqGwsJDeSRFCHELvoDzML37xC6xZs8YtzQkAUlNTsXr1arzxxhtuOR8hRDqoQXmQ0tJSVFRUYPv27W49744dO1BeXo7S0lK3npcQIm7UoDwEwzA4ePAgcnNz4e/v79Zz+/v7Izc3F4cOHQI9USaEsEUNykPodDoMDw8jKSmJl/MnJSXBZDLR0HNCCGvUoDzEu+++i6ysLKfnOc2WXC5HVlYWCgsLeTk/IUR8qEF5AKvViuLiYmRkZPAaR0ZGBs6ePQur1cprHIQQcaAG5QEMBgNCQ0NZra13Pw8++CAOHDjg9OcVCgVCQ0NhNBpnFQchxDNQg/IAer0esbGxfIcBAIiJiYFer+c7DEKICFCD8gANDQ2Ijo7mOwwA1KAIIezRTnMeYGBgABEREQ59hmEYTE5O3vVzm8122zskmUwGLy8v1sdVKBTo6upyKBZCiGeiBuUBJiYmHN71try8fNoJvYcOHcKhQ4em/v/09HS8/fbbrI/r7e0Ni8XiUCyEEM9EDcoD+Pj4ODxybsWKFfjggw9u+9lLL72EzMxMPP3001M/c3R9PavVCl9fX4c+QwjxTNSgPEBISAgGBwcd+kxQUBBWrVp12898fHywYMGCu37uCJPJBKVS6fTnCSGegwZJeACVSoXW1la+wwAAGI1GqNVqvsMghIgANSgPoFarYTAY+A4DADUoQgh71KA8QGxsLPr6+jA0NMRrHCaTCX19fYiJieE1DkKIONA7KA/g7e2NzMxMlJSU4NFHH3X6OJ9//vms4igpKcGmTZscHlFICPFMdAflIXJyclBUVASbzcbL+W02G06fPo3s7Gxezk8IER9qUB5Co9EgODgY1dXVvJzfft6vvvoKBoOB9oUihMyIGpSHkMlkyMvLw7FjxzA2NubWc4+NjaGgoAA/+tGPsHr1ahQXF+PYsWPQ6/W83dERQoRPxtClrEf5yU9+ArPZjBdffNFt5zxy5Aj8/f3x2muvAbi5jFJLSwvOnz+PsbExpKWlYcWKFfRuihByG2pQHmZ4eBjZ2dnIzc1Famqqy89XWVmJY8eOobCw8K5VJxiGwdWrV3H+/Hn09vYiJSUFa9asgZ+fn8vjIoQIHzUoD6TT6fDqq6/ilVdewcqVK112ntraWhw8eBAHDhyARqO577/t7u5GeXk5WltbsWbNGqSkpCAgIMBlsRFChI8alIe6cOEC9u3bh127drnkTqqiogJHjhzB/v37kZyczPpz/f39qKioQENDAxITE5GWlkZLIxHioahBebDS0lL827/9G7RaLXbs2AF/f/9ZH3NsbAzHjx/HV199hfz8/BnvnO5leHgYVVVVuHTpEuLi4pCeno6wsLBZx0cIEQ9qUB5qcnISx48fR0JCAoqKilBeXo4XXngBSUlJkMsdH9xps9lQXV2N3//+9wgPD8eBAwcwd+7cWcc5Pj6OmpoaXLhwAQsXLoRWq8WiRYtmfVxCiPBRg/JQn332GXp7e/Htb38bMpkMpaWlOHToEEwmE7KyspCRkQGFQjHjcUwmE0pKSlBUVASFQoG8vDxcv34dfn5+eOihhziLd2JiArW1tSgvL4dSqYRWq0VMTAxkMhln5yCECAs1KA905coVfPTRR3jhhRcQGBg49XOGYaDT6VBYWIizZ88iNDQUMTExiImJgUKhgLe3N6xWK0wmE4xGI4xGI/r6+rBp0yZkZ2dDo9FAJpNhfHwcR48exZYtWxAfH89p7JOTk2hoaEBZWRnkcjm0Wi2WL1/u1F0fIUTYqEF5GHvzeOihhxAXF3fPf2e1WmE0GqHX66HX6zE4OAiLxQJfX18olUqo1Wqo1WrExMRMO3+pvb0dH374IXbu3HlbE+SKfS5VWVkZRkZGkJ6eTnOpCJEYalAehGEY/PnPf0ZAQAC2bNni8vN9/vnn6O7uxpNPPumyR3H2uVRlZWXo6emhuVSESAg1KA9SV1eHsrIy5ObmwsfHx+Xnm5ycxIkTJ7By5UokJSW5/Hw9PT0oKyubmkuVnJzskrs3Qoh7UIPyEAMDA3jrrbeQk5ODBQsWuO28N27cwIkTJ/Dss88iNDTULeekuVSESAM1KA9gs9nw9ttvY9myZUhLS3P7+S9evIiamho899xzbn1HRHOpCBE3alAe4Ny5c2hvb0dOTg4vw7IZhsEHH3yAuXPnIjMz0+3np7lUhIgTNSiJ6+zsxJ/+9Ce88MILCA4O5i2O0dFRHD16FFu3bkV0dDQvMdjnUlVUVEChUNBcKkIEjhqUhJnNZhQUFCAjIwPLly/nOxwYjUb87W9/w86dOzlZVslZNpsNer2e5lIRInDUoCTsr3/9KwDgG9/4Bs+R/MOZM2cwNDSEb33rW7zfudBcKkKEjRqURDU2NqKkpAQ7d+6Er68v3+FMsVqteOutt5CamurSrT4cZd+XiuZSESIc1KAkaGhoCAUFBfj2t7+NyMhIvsO5S29vL/7whz/gueee42RBWS7RXCpChIMalMQwDIN3330XixcvxoYNG/gO556qqqpQX1+P7du3w8vLi+9w7jIwMICKigro9XqaS0UIT6hBSUxFRQWamprw7LPPCvqlP8MwKCwsRGRkJB544AG+w7mn4eFhXLhwARcvXqS5VIS4GTUoCenp6cE777yD559/HiEhIXyHM6Ph4WEUFBRg27ZtiIqK4juc+zKbzaipqUFVVRUiIiKg1WoFHzMhYkcNSiImJiZw7NgxaLVarFixgu9wWLt8+TLOnDmDnTt3imJQwsTEBOrq6lBeXo7g4GBotVrExsbyPiKRECmiBiURp0+fxtjYGB599FHRfVmWlJRAo9GI6tGZzWab2pdKJpMhPT0dKpVK0I9VCREbalAS0NLSgk8//RQ7d+7EnDlz+A7HYQzDwGazCXKwxEwYhoHBYEBZWRmGh4dpLhUhHKIGJXIjIyMoKCjAo48+iiVLlvAdjktZrVaMjo5idHQUERERfIdzF/u+VF1dXUhJScHatWtF8diSEKGiBiViDMPgT3/6E8LCwrBx40a+w3EZhmFQVVWFM2fOoKenB6Ojo1Cr1UhLS8O6dev4Du8uPT09KC8vh9FoxOrVq5GSkkJzqQhxAjUoEevp6cGnn36KZ555RpSPx9jo6urCr3/9a1gsFsTFxWHZsmUICwvDxYsXceHCBfzoRz9y6/5Wjrh1LpW9oYphdCUhQkENSsQmJibg5eUlyRfzk5OT+Pjjj/H3v/8dK1aswLp166BSqaYWmbVarThy5Ah8fHzw4osv8hzt/Y2MjEztSxUTEwOtViuqASGE8IXe5AoUwzAzjsZzx7btfPnjH/+ImpoaPPTQQ/ja1742tRuvzWYDAHh7e8PLywsymQz2ayyhjl4MDAzExo0bodVqUVNTg3fffZeTuVQMw8BqtWJkZITuzIgkUYMSIKvV6tGjwOrq6lBcXIytW7di69atU3eINptt6j+PjIzg8uXLSEpKEmxjupOfnx+0Wi2Sk5NRV1eHjz/+2Om5VIODgygsLERXVxfGxsYQFRWFp59+GvPmzXPhX0CIe9EjPoH58ssvcf36daSnpwtypJo7fPzxx6isrMRPfvITADcf98nl8qkv8NbWVpw6dQq9vb3YvXu3aL+Ub51LBQBarZbVXKqBgQH87//+L3p7e7Fp0yaEh4ejrKwMg4ODePnllwW3AC8hzvLcy3QBMhqNeP/99/Hwww/zuvst3/r7+xEWFgar1QovL6+pASCjo6Oor69HSUkJuru7sXXrVoSEhLB6HCpEcrkciYmJUKvVU3OpJicnkZiYeN9BLwUFBbBYLNi1a9fURpQLFizAwYMHUVdXJ+i1DQlxBDUoAXn77bexdu1abNiwAYGBgbBarRgbG0NLSwsWLlyI4OBgBAQE8B2my61evRpvvvkmGhoakJiYCACoqamBTqdDZWUl/P398dxzz2HVqlU8R8oNmUyGuLg4xMXF3fYYczqXL1+GwWDAk08+edsuyUuWLEFAQAA6OjoAsHuHSYjQUYMSCIPBALPZjPT09KkmdOrUKdTU1KC/vx82mw0rVqzA5s2boVarJf0FtGrVKqxZswYnTpyAr68vLBYLAMBiseBrX/sasrOzp/6t1PJwv+bEMAwOHz6M1atXY+3atQD+8V6ur68PFosF/v7+kssJ8VzUoAQiODgYZrMZ/v7+kMvlOHfuHL744gtkZWVhyZIlGBgYwJkzZ/C73/0O3/ve9xAfH893yC61Y8cOGI1GVFVVISAgAMHBwVi1atXU8Gz7F7MnfRHrdDrYbDYkJSVNPQK2N7SmpiZMTExg/vz5UyMbPSk3RJqoQQmEr68v5HI5DAYDFi1ahNOnT2Pz5s345je/OfVFo9FocPjwYfz5z3/GP//zP0t6pJ+Pjw+WLVuGZcuW3fZz+5iee91pSPWLmWEYfP7554iMjJy6OLE3aZPJhLq6OgA3B1rci9VqRUtLC4aHhxEdHY358+e7JXZCnCXdbziRUSqViI+PxyeffILo6GjEx8cjNDR06suWYRiEhoZCo9GgoqIC4+PjCAoK4jlq97B/EbNpPvbfz/QuR4xu3LiBpUuXIjg4+La/T6fT4eLFi3jsscfg7e19z7+9s7MTn3zyCQYHBzE6Oor4+Hjk5OR49IAcImzS+l+wyD399NNQKpX42c9+hi+//BJNTU0wm823fTGHh4fDy8sLExMTPEfrPnK5fGrACPCPybp3YhgGk5OT+MMf/oDi4mJ3huhyDMPAYrFM1cGtQ+6/+OILzJ8/H4888shtv7vTkiVL8NJLL2H79u149tlnMTIygl//+te4evWqe/4IQhxEDUogbDYbAgMD8dRTT2HlypUICAjAV199hZKSEvT19WFychLt7e0oLi5GdHS0R8116evrw759+3DkyBEA9368J5PJ4OXlhaamJpw9exafffaZO8N0KblcjgcffBCXL1/G1atXMTo6CpPJhHfeeQcmkwlPPfUUgJt1dK8GxTAM5syZg7i4OKxatQpPPfUUzGYzysvL3fmnEMIaTdQVoO7ublRUVODSpUvo6OjA3LlzERAQgJGREQQGBmLfvn3w9fXlO0y3sVqt+PnPfw6VSoWvf/3r99zCwv5oq7+/H++99x6qq6vx+uuvIzw83M0Ru0Z/fz+OHTuG1tZWREREoLOzEwqFAo8//jiSk5MdOpb9rvzNN99EW1sbXn/9dRdFTYjz6B2UgNi/NBYsWIBHH30Ua9asQXd3N3Q6HQAgJiYGKpXKo5oTcHPdvby8PPj7+993Aqv9zqq/v39qPpBer5dMg5o7dy5+8IMf4KuvvoLRaMSWLVuQkJAAhUIB4B/1c793dXeOfrxy5Qrmz5+P0dFR+Pv7QyaTwWazwWq1YmhoiAZSEF7RHZQATfcFI8WX/lwbGBjARx99hNLSUkRHRyMrKwtxcXGSX0jVXi9msxn9/f33XCLr1rrq6urC6dOnUVVVhcceewybN28GcHNZqYKCArS3t0MulyM4OBjbtm1DTEyM2/4eQuyoQfHk1qtd4O4X23c2KWpQN42OjqK3txdRUVFTd1M2mw3nz5/Hhx9+iMnJSWRmZiI1NRXh4eGSHHJ+LydPnsQXX3yBPXv2TDtPbnR0FM3NzaipqUF1dTWsViuefPJJpKWlTW2o+MEHH+DMmTPQarWIj49HfX09dDodHnnkEWzevBkymcyjckr4RY/4eDA+Pg4/Pz+Mj49jzpw5AG5vQLc2p97eXoSFhVFz+n+VlZUoKirCrl27sHjxYly5cgXvv/8+GhoasHbtWmRkZCA6OnrqPZVU50VNJyQkBKtWrZoaNs4wDIaHh2EwGKaGoo+MjGDx4sXYsmULkpKSEBkZedsxKioq8MADDyAnJwcAkJycjI8//hjnzp1Damqq5O9GibBQg3Kj3t5enDt3DuXl5fDx8UF4eDgSEhLw0EMPTbulRElJCT755BPk5ORMLW3j6TZu3IiPPvoIn376KZRKJUpKShAaGornnnsOGo3mri9QT2lOALB169bbtmo5ffo0Pv/8c4yPj0OpVGLdunVYs2YNoqKibttL7NZ5ZosWLcLIyMjU7+wTpv/+97+joaHhvhOBCeEaNSg3Onz4MGQyGVJSUjA5OYkrV67gk08+wZkzZ5CTk4Pk5OTbJqQGBATAZDJN3WWRm55//nn85je/gbe3NzZu3Ij169djwYIFU1/Md941jY6Ooq2tDf7+/oiOjuYpavfw9vaemjN1/fp19PX1YenSpcjNzb3tbunWR8tyuRyTk5NoampCREQE6uvr0draiujoaAwNDaG7uxsAsHTpUl7+JuK56B2Um3z55Zf4+OOPsWfPnqndYU0mE/R6Pc6fP4+GhgakpqYiJyfnthUibn0MSP7h8OHDMBqNePnll2fclXZ8fBwHDx6E2WzGww8/jNTUVDdFyb+mpiacPHkSQ0ND2LhxIzZu3HjXXea1a9dQWFiI9vZ2zJ07F35+fmhuboZKpcLk5CTa2tqg1WrxzDPPeNQdKeEfNSg3+eijj6DT6fDqq69izpw5YBhm6lFed3c3zp8/j7KyMqhUKuTk5MDf35/niIVtaGgI//Iv/4ItW7Zg69attz2yupX9bqqzsxMVFRUoLS3Fv//7v08NzfYU58+fx/vvvw+ZTIbHHnsMWq0W3t7esFgs+Oijj3D69Gns2bMHixcvBgA0NDTgnXfeQVpaGtauXSv5O08iTPTm3U2USiU6OjpgtVqnHqvYl+xZsGABtm7digceeADl5eWora3lOVrhCw4OxsMPP4y+vr77XtXbfxcZGTm1Mvx7773nrjAFY926dcjPz8fmzZtRW1s7lRcfH5+pASbx8fHw8/ODn58fVqxYgdjYWBgMBixduhQMw4CuZYm7UYNyk5SUFISEhOC3v/0trly5AgBT75tsNht8fHywdetWqFQqVFVV8RytODz22GN48cUXWa/qHhgYCLVaja6uLphMJhdHJ0wPP/wwvve978HLywsMw2BiYgKjo6MYGxub2neLYRh4eXlh8eLF6O3txeTkJA0vJ7ygBuUmAQEByMnJQX9/P/70pz+hpqYGZrP5trsp+yiqsbExmM1mvkMWjWvXrsFgMKCnpwfj4+N3/f7WAQF9fX30Xu//yWQy+Pr64hvf+Ab6+vqg0+mmFiEeHx/HF198gaioqLtq0d7YCHE1GsXnRitWrMDu3bvx61//GidOnIBWq8XKlSuhUqkgl8tx48YNGAwGREZG3nO9OXI3g8GA48ePw9/fH35+flPblYSHhyMiIgKBgYGQyWSoqqpCVVUVNBqNpPfSclRSUhLa2tpw9OhRrF69Gl5eXmhubsbY2BjS09OnJvHa2Ww2/PGPf0RgYCC0Wq1klpIiwkODJHjAMAxKS0vxzjvvQKlUIigoCGFhYejo6IDFYsFrr73mcevtzYbVasXevXsRFRWFuLg4GAwGdHd3Y3R0FDabDUFBQRgfH4eXlxcUCgW2b9+OhIQEvsMWHIPBgDNnzmB0dBRyuRxZWVlQqVTTPtozm824ePEiqqqqEB4eDq1Wi6ioKHoMSDhFDYpH9scoBoMB169fR3JyMhITE6dGUhH2zp07h/fffx//+q//itDQUExMTKC7uxs9PT3o6uqa2rF4ujsCcjuz2cz6Dt5qtaKurg4VFRUICAiAVqtFXFwcNSrCCWpQRDJ++tOfQqlU4sUXX7zvHSita8g9m82GxsZGlJWVwWazQavVQq1WU57JrFCDchP7Fgb06M51urq68Prrr+O73/3u1NJQbLagINxhGAZGoxHl5eUYHBxEWloaVq5cec95aoTcDzUoN2AYBh9++CGWLVuGxMREvsORtJMnT6KyshL//d//jYCAAL7D8WgdHR0oLy9HZ2cnkpOTsXbtWho9SRxCDcoNamtrUVlZieeff55Gj7lBTU0NLa4rIL29vaioqEBzczNWr16NlJSU25bzIuReqEG5WH9/P44fP45nnnkGYWFhfIdDCG8GBwdRWVkJnU4HlUqFtLQ0zJ07l++wiIBRg3Ihm82GkydPQq1WIyUlhe9wCG6OOvPy8qL3UTwaHR1FVVUVLl68iOjoaKSnp2PBggV8h0UEiBqUC3355Zfo7OzEU089RV+IAsAwDK5cuYLW1lZs3LiR73A8Hs2lIjOhBuUiV69exalTp/DCCy/Q83YBGRsbw9GjR/H1r38dMTExfIdDcPOuVqfToby8HP7+/li3bh3NpSIAqEG5hNlsxtGjR7F582ZasUCAWltb8Ze//AU7d+6kkX4CYrPZ0NTUhLKyMkxOTkKr1UKlUsHLy4vv0AhPqEG5wF/+8hd4eXnhkUce4TsUcg9nz57FwMAAHn/8cbpSFxiGYdDa2oqysjIMDAwgLS0Nq1atorlUHogaFMf0ej2+/PJL5Obm0qRcAbNarTh+/DiSkpKwevVqvsMh93DrXKqkpCQkJSXRXCoPQg2KQyaTCceOHcOTTz6JhQsX8h0OmUFfXx/efvtt7NixA/PmzeM7HHIffX19KC8vR3NzM1atWoXU1FR6t+sBqEFxxGaz4d1330V0dDTWrVvHdziEpQsXLkCn02H79u30rkMEbp1LtXz5cqSnp9NcKgmjBsWR8vJytLS04Omnn6YFMkWEYRj88Y9/xIIFC/Dggw/yHQ5haXR0FBcuXEBNTQ2WLl0KrVZLc6kkiBoUB7q6uvDee+8hNzcXCoWC73CIg0ZGRnD06FF861vfoq1ORMZsNuPSpUuorKxEWFgYtFotFi9eTANfJIIa1CxNTEygoKAAGzZsoIVgRay5uRmnT5/Gd77zHdrNWITunEul1WoRHx9PjUrkqEHN0qeffgqLxYJvfvObfIdCZunTTz+F2WzGo48+yncoxEl3zqVKT0+HWq2m94siRQ1qFpqbm1FUVISdO3fSVbcE0N2wdNjnUpWXl6O/v5/mUokUNSgnDQ8Po6CgANu2bUNUVBTf4RCO2N8nPv/881AqlXyHQzjQ2dmJ8vJyXL16FcnJyTSXSkSoQTmBYRi89957iIiIoJFfEmSfb/PMM8/QiEwJuXMuVUpKCoKDg/kOi9wHNSgn0NwZaWMYBu+88w7NaZMok8mEioqKqblUaWlpNFFboKhBOai3txd/+MMfaPUBiaNVQaSP5lIJn8c1KKvVCoPBAL1ej4aGBgwMDGBiYgI+Pj4ICQmBSqWCWq1GbGzsXduzW61WvPXWW0hOTqb12zyAXq/HF198gRdeeGHadRVnU0tEOCwWCy5evMj7XCqqp7t5RINiGAZ1dXUoLCxEcXExQkNDERsbi+joaCiVSnh7e8NqtWJwcBCtra0wGAzo6+vDpk2bkJ2dDY1GA5lMhrNnz2JwcBDbtm2j+RUe4s6V6bmqJSI89rlUFRUVmDNnjlvmUlE93Z/kG1RpaSkOHjyI4eFhZGVlISMjg9WL0aGhIRQXF6OoqAgKhQJPPfUUOjs78Z3vfAf+/v5uiJwIgdlsRkFBATIzM9HX18dJLeXl5WH9+vVuiJ44wz6Xqry8HBMTE9BqtS6ZS8XVd5OU60myDWpoaAhvvPEGKioqkJubi6SkJKdGZNlsNlRXV+N3v/sdVq9ejR//+Me0irKHaWxsxOuvv47+/n688MILs66lgoICpKenY/fu3VRLAsYwDNra2lBWVsbpXCquv5ukXE+SbFA6nQ579+7FmjVrsH37dk7ueMbGxnDy5ElcvHgR+fn50Gg0HERKhM5eS6tWrcJzzz1HteShuJpLRd9NjpFcg7pw4QL27duHXbt2ITU1lfPjV1ZW4vDhw9i/fz+Sk5M5Pz4RDqolcqe+vj5UVFTg8uXLWLlyJVJTU1nPpaJ6cpykGpROp8Orr76KV155BStXrnTZeWpra3Hw4EEcOHBAUlcr5B+olsj9mEwmVFZWoq6uDsuWLUN6evp9p51QPTlHMg1qeHgY2dnZyM3NdcnVyZ0qKytx7NgxFBYWSu65r6ejWiJsjY6Oorq6GtXV1ViyZAm0Wi0iIiJu+zdUT86TTIP6yU9+ArPZjBdffNFt5zxy5Aj8/f3x2muvue2cxPWoloijLBbL1L5U8+fPh1arxZIlSyCTyaieZkESC42VlpaioqIC27dvd+t5d+zYgfLycpSWlrr1vMR1qJaIM3x9fZGamoqXXnoJarUan376KY4fP4733nuP6mkWRH8HxTAMnn32WTz++ONISUlx+/mrqqpw6tQpnDx5UtIT5jwB1RLhin0u1Z49e/D8889TPTlJ9HdQOp0Ow8PDSEpK4uX8SUlJMJlM0Ol0vJyfcIdqiXBFLpfDarVCJpNRPc2C6BvUu+++i6ysLN62RZDL5cjKykJhYSEv5yfcoVoiXKJ6mj1RNyir1Yri4mJkZGTwGkdGRgbOnj0Lq9XKaxzEeVRLhEtUT9wQdYMyGAwIDQ11aNOx0dFRbNmyBY8//jgmJiamfv7FF18gPj4ex48fdzgOhUKB0NBQGI1Ghz9LhMGZWgK4ryeqJWmgeuKGqBuUXq9HbGysQ58JCAjAr371KzQ0NOCXv/wlgJuzw/fs2YPNmzfjueeecyqWmJgY6PV6pz5L+OdMLQGuqSeqJfGjeuKGqBtUQ0MDoqOjHf6cRqPB3r17ceTIEZw7dw579uyBl5cXfvrTnzodi5iLgDhfSwD39US1JH5UT9wQ9a5XAwMDd83aZmvnzp348ssv8eKLL2JiYgLHjh2b1Q65CoUCXV1dTn+e8Gs2tQRwW09US+JH9cQNUd9BTUxMOL2zpEwmw7Zt22CxWKBWq7Fhw4ZZxeLt7Q2LxTKrYxD+zKaWAG7riWpJ/KieuCHqBuXj4+P06JTe3l7853/+JzQaDfR6PY4ePTqrWKxW67TbghNxmE0tAdzWE9WS+FE9cUPUDSokJASDg4MOf45hGOzduxe+vr44fvw4du7cif3796OhocHpWEwmE5RKpdOfJ/xytpYA7uuJakn8qJ64IeoGpVKp0Nra6vDn3nzzTZw7dw5vvPEGlEol9u7di4SEBPzgBz/A+Pi4U7EYjUao1WqnPkv452wtAdzXE9WS+FE9cUPUDUqtVsNgMDj0mbq6OvziF7/AP/3TPyE9PR3AzYUef/WrX+Hq1av4r//6L6diEXMREOdqCXBNPVEtiR/VEzdEPYovNjYWfX19GBoaYj0hbsWKFdMOuYyNjUVdXZ1TcZhMJvT19SEmJsapzxP+OVNLAPf1RLUkDVRP3BD1HZS3tzcyMzNRUlLCaxwlJSXYtGnTrEbtEH5RLREuUT1xQ9QNCgBycnJQVFQEm83Gy/ltNhuKioqQnZ3Ny/kJd6iWCJeonmZP9A1Ko9EgODgY1dXVvJy/uroaSqUSGo2Gl/MT7lAtES5RPc2e6BuUTCZDXl4ejh07hrGxMbeee2xsDAUFBcjLyxPthmDkH6iWCJeonmZP9A0KANavX4+0tDScPHnSrec9ceIEtFot1q1b59bzEtehWiJconqaHUk0KADYvXs3Ll68iMrKSrecr7KyEpcuXcIPf/hDt5yPuA/VEuES1ZPzJNOggoKCkJ+fjyNHjqC2ttal56qtrcWRI0eQn5+PoKAgl56LuB/VEuES1ZPzZAzDMHwHwaULFy5g37592LVrF1JTUzk/fkVFBY4cOYL9+/cjOTmZ8+MT4aBaIlyienKc5BoUAOh0Ouzduxdr1qzB9u3b4e/vP+tjjo2N4cSJE7h06RLy8/NFPTKGsEe1RLhE9eQYyTziu5VGo0FhYSH8/Pywd+9eVFVVOT0XwWazoaqqCnv37oW/vz8KCwslVQDk/qiWCJeonhwjyTuoW5WWluLQoUO4ceMGHnnkEWRkZEChUMz4OZPJhJKSEhQVFUGhUCAvLw/r1693Q8REqOy1NDg4iKysLGRmZrKupeLiYpw5c4ZqiUyx19PAwAAeeughh76bPKWeJN+gAODq1asoKCiA2WxGcXExQkNDERMTg5iYGCgUCnh7e8NqtcJkMsFoNMJoNKK3txerVq3Cyy+/DI1GI+q5BIQ7DMPgnXfeQXFxMerr61nVUl9fH+Li4rB9+3Zs2bKFaolMYRgGv/rVr9DW1oaqqirW9bR48WLs27cPa9askXQ9iXOBJgfp9XpkZmZiw4YNsFqtMBqN0Ov10Ov16OrqgsViga+vL5RKJTIzM/H9738fQUFBOHXqFDUncpeRkRG89tpriIiIYFVLMTExqKmpQU9PD9USuc3IyAhkMhny8/Mhk8lY19OHH34Ib29vydeT5BvU5OQkGhoasGPHDgA3F3FMSEhAQkICHnvssXt+jmEY+Pr64urVq1i8eLG7wiUC19vbi4mJCSxatAgymYxVLQE3t184d+4cJiYm4OPj46ZoidDp9XokJCRM1QTbekpMTMRXX32FlStXuiNM3khykMStWltboVQqMXfuXIc+J5PJkJiYiPr6ehdFRsRIp9MhMTHR4SvXoKAgLFy4EM3NzS6KjIhRfX09EhMTHf5cfHw8rl27huHhYRdEJRySb1B6vd7pkS2JiYlobGzE5OQkx1ERMWIYBnq93qkvFAB0wUNuc/36dQwNDWHp0qUOf9bHxwcJCQlObwMvFpJuUBaLBc3NzVCpVE59XqlUYt68eTAajRxHRsSovb0dc+bMQVhYmFOfX7ZsGa5cueL2hUOJMNXX10OlUkEud+5rWK1WS/6CR9INqrm5GZGRkQgMDHT6GBqNBjqdjsOoiFjV19fPap6Jn58fYmJi0NjYyGFURIwYhpl1PUVHR8NkMuHGjRscRiYskm5Qzj7fvdXy5cthMBhgNps5ioqIkdVqRWNjI9Rq9ayOo9FoJH/VS2Z27do1yGQyREREOH0MuVwOlUol6XqSbIMaHR1Fe3s7EhISZnWcgIAALF68GJcvX+YoMiJGRqMRYWFhrCZS3k9MTAx6e3thMpk4ioyIkf3iebbDxO3vNaU6nVWyDaqxsRGxsbHw8/Ob9bHo5Taxj96bLW9vbyxfvhx6vZ6DqIgY2Ww2NDQ0cFJPCxcuBMMw6Orq4iAy4ZFsg9LpdJytSxUfH4/Ozk6MjIxwcjwiLmazGUajEcuXL+fkeImJifRe04O1trZCoVBg3rx5sz6WfTqMVOtJkg1qcHAQN27cQExMDCfH8/X1RXx8vOSHdJLpNTU1YenSpZysPA0Aixcvxvj4OHp7ezk5HhGX2UxVmE5iYiIaGhqcXnRWyCTZoOrr67F8+XJ4eXlxdky1Wi3ZqxRyf/X19bMeHHErmUzmEUOEyd0mJiZw+fJlp6e+TGf+/PkIDg5GW1sbZ8cUCsk1KPvwTS6vUICbL7cHBwfR39/P6XGJsA0PD+PatWuIj4/n9Lj20XxSfblNptfc3IyFCxdyvtutVN+TS65B9fb2wmKxICoqitPjesKQTnK3O9dK40pYWBh8fHzQ0dHB6XGJsLni4hm4+YTn8uXLmJiY4PzYfJJcg7I/jnHFKr9SH9JJ7sb1+wI7mUxGc6I8zNjYGK5cuYJly5ZxfuygoCBERESgpaWF82PzSVINiovZ2fcTGRkJm82G7u5ulxyfCMuNGzdgMpmcWiuNDfvLbVrr0TM0NjYiJiaGk6kv05HiYz5JNairV6/Oaq20mdDLbc8y27XSZkJrPXoWV148AzdXvWlra8P4+LjLzuFukmpQXE2mvJ/ExETo9XpJDukk/+Dqu3E7KV71kruZTCb09vYiNjbWZeewr/UopekwkmlQk5OTaGpq4nQ48HRCQ0MRGBiIK1euuPQ8hF/2mfmzWSuNDZVKBYPBAIvF4tLzEH65YurLdKR2wSOZBmUwGDB//nwolUqXn0tqRUDu5uzGhI4KCAhAVFQUrfUoca4avXen2NhYSa31KJkG5Y7HMXZqtRpNTU2wWq1uOR9xLy7XSmNDykvVkJtTX8bHx7F48WKXn8vb2xvLli2TzFqPkmhQZrMZBoOBs7XSZhIcHIwFCxZIbkgnuamtrQ3BwcGcrJXGBq31KG2unPoyHSk94ZFEg2pqasKSJUs4WyuNDSkVAbmdux7H2Pn6+iIuLk5SL7fJTe4abHOrJUuWYHR0VBJrPUqiQbn7CwW4OaSztbVVUkM6yT/WSnP1YJs70QWPNHV0dMDX19dlU1+mY1/hXAqP+UTfoFy1VtpM5syZg+joaNq+W2JaWlpcslbaTKKjozEwMICBgQG3npe4lrsG29xJKqveiL5BNTQ0ID4+nvO10tigq17p4eNuHAC8vLxorUeJmZycRGNjIy/1FB4eDm9vb9Gv9Sj6BsXXFwoAxMXFoaenB0NDQ7ycn3BrbGwMbW1tLlkrjQ37aD6xX/WSm4xGI+bNm+eWqS93sj/mE/sFj6gblH2ttOjoaF7O7+3tjYSEBEk86yWuXyttJpGRkZicnKS1HiXC3YMj7pSYmIjGxkZRr/Uo6gbl6rXS2JDCVQq5ic+7cUA6V73k5tSXlpYWt019mU5ISAhCQkLQ2trKWwyzJdoG5aqNCR21ZMkSjIyMoK+vj9c4yOy4Y600NmitR2m4fPkyFi9ejICAAF7j0Gg0op4ELtoG1dXVBYZhsHDhQl7jkMvltMK5BOj1eixfvhze3t68xhEaGoqAgAC0t7fzGgeZHSFcPAM313psaWkR7VqPom1Q9gJw9/DN6diveunltngJ5QsFoMfGYjcyMoLOzk63T32ZjtjXehRlg7LZbC7b6dQZCxYsgFwuR2dnJ9+hECf09vZibGzMLWulsUFrPYpbQ0MD4uLi4Ovry3coAMR9wSPKBnXlyhUEBwdj/vz5fIcCgF5ui52710qbiUKhQHh4OK31KFJCuhsHgISEBHR0dGB0dJTvUBwmygbljo0JHUXbd4sTH2ulsUEXPOLU39+PgYEB3qa+TEfMaz2KrkHZ10pTqVR8h3KbuXPnQqlUinpIpyfq6OiAj4+PW9dKY4PWehQn+9QXV29M6CixbukiugbV0tKCiIgIBAcH8x3KXTQaDU3aFRn73ZNQHu/ZzZkzB0uXLkVTUxPfoRCWhDL1ZTpiXetRdA1KqAUA3BzS2dzcLNohnZ5mcnISDQ0Nbl+5nC16zCcu3d3dmJycRGRkJN+h3MXLywvLly8XXT2JqkGNj4+jra2N19nZ9xMYGIjIyEg0NzfzHQphobW1FfPmzUNISAjfoUwrPj4e3d3dtNajSAhp6st0xLjCuagaVENDA6Kjo3lbK40NuuoVDyEOtrkVrfUoHkKb+jKdRYsWwWq1oqenh+9QWBNVg9Lr9YIbbXWnhIQEtLe3i3JIpyexWCwwGAyCG2xzJ7rgEYf29nYEBgYiNDSU71DuSYzTYUTToEwmE3p6enhfK20mfn5+iIuLo40MBe7y5cuIiorifa20mdjXerx+/TrfoZD7EPrduJ29QYllrUfRNCi9Xo9ly5bxvlYaG2Id0ulJhDzY5lZyuZw2MhQ4q9WKy5cvC3awza3EttajaBqUWL5QACAmJgY3btzA4OAg36GQaYyMjKCjo0MQa6WxIcaX256kpaUF4eHhgpz6Mh0xPeYTRYPq6+vD6OioYNZKm4lYh3R6CqGtlTaTiIgIWutRwMR08QyIa61HUTQo+1ppfG5M6Ci66hUusX2h2F9u02g+4RkfH0dra6tgp75MR6FQICwsDAaDge9QZiT4b3yhrpU2k6ioKFgsFvT29vIdCrmFfTa9kNZKY4M2MhSmxsZGREdHY86cOXyH4hCxvCcXfIPq6OiAt7c3wsPD+Q7FIWIc0ukJhLpW2kxorUdhEvrcp3tRqVRobW2F2WzmO5T7EnyDEvrs7Puhx3zCwjCMaIYDT4cueIRlaGgI3d3diIuL4zsUh9nXehT6dBhBN6jJyUk0NjaK9gslLCwMc+bMEc2QTqnr6ekR7FppbKjVajQ3N2NiYoLvUAhu3j0lJCSIYurLdMRwwSPoBtXa2oqQkBDBrpXGhhiKwFPY757EeDcO3FzrceHChaLdvltqxDbY5k5xcXHo6urC8PAw36Hck6AblBgHR9wpMTERTU1NtJEhz8SwVhobGo2GLngE4Pr16xgZGcGSJUv4DsVpPj4+gl/rUbANymKxoKWlRVTDN6ejUCgwf/58UQzplLL29nYEBAQIeq00NuxrPY6NjfEdikfT6XRQqVSimvoyHaE/4RFsdpubm7Fo0SIEBgbyHcqs0VUv/8T+OMbOz88PsbGxoty+WyoYhhHFwtVsLF26FENDQ4Jd61GwDUrMo63utHz5chgMBsEP6ZQqq9WKpqYmUayVxobQr3qlrrOzE3K5HAsWLOA7lFmTy+VQq9WCrSdBNqjR0VF0dHQgISGB71A44e/vjyVLltD23Tyxr5WmUCj4DoUTsbGxuH79Oq31yBMxT32ZjpCnwwiyQTU0NCA2NlY0a6WxQVe9/JHK4z07Ly8vLFu2TNAvt6VqcnISDQ0NkqqniIgIyGQyXLt2je9Q7iLIBiWF0Xt3io+Px7Vr1wQ9pFOKzGaz6NZKY0Oj0YhiqRqpaWtrg1KpxNy5c/kOhTNCXvVGcA1qYGAAN27cEN1aaTOxD+mkl9vuJda10mYSFRUFs9ksqu27pUCKF8/AzSc8DQ0NglvrUXANSqxrpbEh5JeRUmVfCV9qhHzVK1UWiwXNzc1QqVR8h8K5efPmQaFQCG6tR0E1KPvK5VJ6vnur6OhomEwm3Lhxg+9QPIJ9rTSxbEzoKCG/3Jai5uZmREZGSmLqy3SEuKWLoBpUT08PrFYrFi1axHcoLkHbd7uX2NdKm0l4eDj8/Pxw9epVvkPxCFK+eAZurnB++fJlQa31KKgGZX8cI5Xhm9Ohq173kfoXCkCDJdxldHQU7e3tkpn6Mp2goCAsXLgQzc3NfIcyRTANSuqP9+wWLlwIhmHQ1dXFdyiSJoW10tiwb99Naz26VmNjI2JjY+Hn58d3KC4ltPeagmlQ9rXSwsLC+A7FpejltnvYB9uIfa20mSiVSlrr0Q2kOnrvTsuWLRPUWo+C+V+vlJY2mglt3+1a9rtxT/hCAYR31Ss1g4ODuH79OmJiYvgOxeX8/PwQExMjmI0MBdGgpLZW2kzmz5+P4OBgtLW18R2KJF27dk0ya6WxQWs9ulZ9fT2WL18uyakv00lMTBTMe01BNCiDwYCwsDDJrJXGBl31uo7U1kqbSUBAABYvXkwbGbqIJ7wbv5V9rUeTycR3KHDp+Fur1QqDwQC9Xo+GhgYMDAxgYmICPj4+CAkJgUqlglqtRm1trUcVAHDz5fa5c+cwMTEBmUzGKk+xsbGSHTI9E7a1FB0dDb1ejx07dvAdslslJiaitrYWKpWKaokFtvUUFBQEi8WCqKgovkN2G/taj/X19UhJSeG1nmQMx+OdGYZBXV0dCgsLUVxcjNDQUMTGxiI6OhpKpRLe3t6wWq0YHBxEa2srDAYDuru7sWnTJjzzzDPQaDQeceXLMAx++ctf4sqVK6iqqmKVp76+PmzatAnZ2dkekSdnaqm3txdxcXHYt2+fR+QIuJmnixcv4mc/+xmuXr1KtXQPztbTypUrkZeX51F5KioqwltvvQWj0chrPXHaoEpLS3Hw4EEMDw8jKysLGRkZCA4OnvFzQ0NDKC4uRlFRERQKBfLy8rB+/XquwhIce55MJhO2bNmCzMxMytMdqJbYuTVPmzdvplq6B6ondoRWT5w0qKGhIbzxxhuoqKhAbm4ukpKSnBrea7PZUF1djYKCAqSnp2P37t0ICgqabXiCQXmaGeWIHcoTO5QndoSap1kPktDpdMjJyYHZbMb+/fuRkpLi9NwTuVyOlJQU5Ofnw2w2Izs7WzCjSWaL8jQzyhE7lCd2KE/sCDlPs7qDunDhAvbt24ddu3YhNTXV6SDupbKyEocPH8b+/fuRnJzM+fHdhfI0M8oRO5QndihP7Ag9T043KJ1Oh1dffRWvvPIKVq5c6cwhWKmtrcXBgwdx4MABUU68pDzNjHLEDuWJHcoTO2LIk1MNanh4GNnZ2cjNzXVJ171TZWUljh07hsLCQlE996U8zYxyxA7liR3KEztiyZNTDxp/8YtfYM2aNW75wwAgNTUVq1evxhtvvOGW83GF8jQzyhE7lCd2KE/siCVPDjeo0tJSVFRUYPv27Y5+dFZ27NiB8vJylJaWuvW8zqI8zYxyxA7liR3KEztiypNDDYphGBw8eBC5ubnw9/d3OMDZ8Pf3R25uLg4dOiT4vZQoTzOjHLFDeWKH8sSO2PLkUIPS6XQYHh5GUlLSbT9/6aWXsHbt2nsuVjk8PIwVK1Zg7969uHbtGv7jP/4DTz75JDQaDeLi4ljvCJqUlASTyST44Z1c5Olvf/sbvv/97+OBBx5AYmIisrKykJ+fj+Hh4RnPL4Y83StHgHvqSQw5ArjLk7P15Gl5onoSVj051KDeffddZGVl3TVG/oknnoDJZEJxcfG0n/vkk08wNjaGJ554Am1tbfjrX/8KhULh8PNPuVyOrKwsFBYWOvQ5d+MiT7/73e/g5eWFPXv24OjRo9i+fTtOnjyJ3NzcGbfpEEOe7pUjwD31JIYcAdzlydl68rQ8UT0Jq55YNyir1Yri4mJkZGTc9bvMzEzMnTsXH3zwwbSfff/99xEZGQmtVou0tDRUVFTg97//Pb7+9a+zPf2UjIwMnD17Flar1eHPugNXeTp8+DAOHjyIb33rW0hPT8fOnTvx2muv4eLFiygrK5sxDiHn6X45AtxXT0LOEcBtnmZTT56UJ6onYdUT6wZlMBgQGho67bpMvr6++OY3v4nPPvsM/f39t/2us7MTFRUV2LZtG2Qy2ax3OFUoFAgNDYXRaJzVcVyFqzzNnz//rs+vWrUKANDd3T1jHELO0/1yBLivnoScI4DbPM2mnjwpT1RPwqon1v9t6PV6xMbG3vP3TzzxBCYmJvDxxx/f9vNTp06BYRg88cQTbE81o5iYGOj1es6OxyVX5qmiogIAEBcXxyoWoeZpphwB7qsnoeYIcH2eHKknT86TIzw5T66oJ9YNqqGhAdHR0ff8/apVq5CQkHDX7eGpU6ewdu1aTrdLFnIRuCpPXV1d+OUvf4kNGzZMXanMRKh5milHgPvqSag5AlybJ0fryVPz5ChPzZOr6ol1gxoYGIBSqbzvv3n88cdx6dKlqVu3S5cuoaWlhdO7J+DmLeLg4CCnx+SKK/I0MjKCl156Cd7e3vjZz37GOhah5olNjgD31JNQcwS4Lk/O1JMn5skZnpgnV9YT6wY1MTEx446J27Ztg1wux/vvvw8A+OCDD+Dr64utW7eyPQ0r3t7esFgsnB6TK1znaXx8HN/97nfR3t6OgoICLFy4kHUsQs0TmxwB7qknoeYIcE2enK0nT8uTszwtT66uJ9YNysfHZ8ZRFwsWLMCGDRvw4YcfwmKx4C9/+Qs2b97Mqms7wmq1wtfXl9NjcoXLPE1MTODll19GXV0d3nzzTSxfvtyhWISaJzY5AtxTT0LNEcB9nmZTT56Up9nwpDy5o55YN6iQkBBWt2RPPPEEOjo68POf/xw3btzg/PEeAJhMJs6bHle4ypPNZsMPf/hDnD9/Hr/97W+xdu1ah2MRap7Y5ghwfT0JNUcAt3mabT15Sp5my1Py5K56mvl+7/+pVCp89tlnM/67hx56CEFBQfj973+P+fPn48EHH7zr3/ztb38DANTV1QEAPvvsM8ybNw/z5s1Denr6jOcwGo3IzMxkG7pbcZWnH//4x/jrX/+K73//+wgICEBNTc3U7yIiIljdSgs1T2xzBLi+noSaI4DbPM22njwlTwDVEyCcemJ9B6VWq2EwGGb8d3PmzME3vvENMAyDxx57bNpnnnl5ecjLy8Pbb78NAHjttdeQl5eHAwcOsIrFaDRCrVazDd2tuMqTvZB+85vf4Mknn7zt/9jOwhZqntjmCHB9PQk1RwC3eZptPXlKngCqJ0A49cR6Pyir1YqMjAz8z//8zz0nermDyWRCXl4eSkpKWL3wczfK08woR+xQntihPLEjxjyxvoPy9vZGZmYmSkpKZhvfrJSUlGDTpk2CLACA8sQG5YgdyhM7lCd2xJgnh9b1yMnJQVFR0YyLlbqKzWZDUVERsrOzeTk/W5SnmVGO2KE8sUN5YkdseXKoQWk0GgQHB6O6utqp4GaruroaSqXS4X3t3Y3yNDPKETuUJ3YoT+yILU8ONSiZTIa8vDwcO3YMY2NjTgXorLGxMRQUFCAvLw8ymcyt53YU5WlmlCN2KE/sUJ7YEVueHF66d/369UhLS8PJkycdDnA2Tpw4Aa1Wi3Xr1rn1vM6iPM2McsQO5YkdyhM7YsqTU2vL7969GxcvXkRlZaUzH3dYZWUlLl26hB/+8IduOR9XKE8zoxyxQ3lih/LEjljy5FSDCgoKQn5+Po4cOYLa2lpnDsFabW0tjhw5gvz8fAQFBbn0XFyjPM2McsQO5YkdyhM7YskT63lQ07lw4QL27duHXbt2Obx9OxsVFRU4cuQI9u/fj+TkZM6P7y6Up5lRjtihPLFDeWJH6HmaVYMCAJ1Oh71792LNmjXYvn07/P39Z3M4ADdfpp04cQKXLl1Cfn6+4EfGsEF5mhnliB3KEzuUJ3aEnKfZ7b+Om8MWCwsL4efnh71796KqqsrpMfY2mw1VVVXYu3cv/P39UVhYKIkCAChPbFCO2KE8sUN5YkfIeZr1HdStSktLcejQIZhMJmRlZSEjIwMKhWLGz5lMJpSUlKCoqAgKhQJ5eXlYv349V2EJDuVpZpQjdihP7FCe2BFanjhtUADAMAx0Oh0KCwtx9uxZhIaGIiYmBjExMVAoFPD29obVaoXJZILRaITRaERfXx82bdqE7OxsaDQawc8l4ALlaWaUI3YoT+xQntgRUp44b1C3slqtMBqN0Ov10Ov1GBwchMViga+vL5RKJdRqNdRqNWJiYgS7fpU7UJ5mRjlih/LEDuWJHb7z5NIGRQghhDhr1oMkCCGEEFegBkUIIUSQqEERQggRJGpQhBBCBIkaFCGEEEGiBkUIIUSQqEERQggRJGpQhBBCBIkaFCGEEEGiBkUIIUSQqEERQggRJGpQhBBCBIkaFCGEEEGiBkUIIUSQqEERQggRJGpQhBBCBOn/AKGKQEWkNQ3rAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from spn.io.Graphics import plot_spn\n",
    "plot_spn(spn_marg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c9476626-7ee1-47d3-8641-d3ba7d678a53",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABOk0lEQVR4nO3deVRT57o/8G8ggCAkICjiDMoQI4IgguDApHWoE7XYXmwptx49bemyR8V11h16em7vWb2VauvB01pRC3WOdeqp2ioKVsUCgiBEoApBBWsVBUIgQEL27w9/5DigsEOmHZ7PWv4hsPf75rvfnSf7zR54DMMwIIQQQsyMlak7QAghhHSHChQhhBCzRAWKEEKIWaICRQghxCxRgSKEEGKW+KbuACHEuBiGQVlZGSQSCbKzs+Hm5gYvLy+MGTMGQqEQfD4farUaTU1NqKmpQXV1Nerr6xEdHY34+HiIxWLweDxTvwzSD/DoNHNC+o/c3FykpaVBoVAgNjYWkZGRcHJy6nG55uZmZGdnIysrCwKBAMnJyQgPDzdCj0l/RgWKkH6gubkZmzZtQn5+PhITExEUFAQrK/Yz/BqNBkVFRcjIyEBoaCjWrl0LR0dHA/SYEPoOihCLJ5VKsWzZMrS3t2PDhg2YPHmyTsUJAKysrDB58mSkpqaivb0d8fHxkEqleu4xIY/QERQhFqywsBDr16/HypUrERISovf1FxQUYNu2bdiwYQOCg4P1vn7Sv1GBIsRCSaVSrF69Gu+//z78/f0N1k5paSnS0tKwefNmiMVig7VD+h8qUIRYIIVCgfj4eCQmJhrkyOlpBQUFyMzMhEQioe+kiN7Qd1CEWKCNGzciMDDQKMUJAEJCQhAQEIBNmzYZpT3SP1CBIsTC5ObmIj8/HwkJCUZtd/ny5cjLy0Nubq5R2yWWiwoUIRaEYRikpaUhMTER9vb2Rm3b3t4eiYmJ2LJlC+ibA6IPVKAIsSBSqRQKhQJBQUEmaT8oKAhyuZxOPSd6QQWKEAty4MABxMbG6nydU19ZWVkhNjYWEonEJO0Ty0IFihALoVarkZ2djcjIyD6tZ8aMGdi8ebPOy0dGRuLs2bNQq9V96gchVKAIsRDV1dVwc3Pr1b31DEkgEMDNzQ0ymcyk/SDcRwWKEAtRXl4OLy8vU3cDAODp6Yny8nJTd4NwHBUoQixERUUFxowZw2oZhmGgVquf+Ac8uins4z/r7OxktV4qUEQf6HlQhFiIxsZGDB06lNUyeXl53V4vtWXLFmzZskX7/9DQUOzdu7fX6xUIBLh79y6rvhDyNCpQhFgIlUoFPp/dLj1hwgQcOXLkiZ+tWrUKUVFReO2117Q/Y3v7Ij6fj46ODlbLEPI0KlCEWAgbGxvWZ845Ojpi4sSJz6zH3d39mZ+zoVarYWtrq/PyhAD0HRQhFsPZ2RlNTU2m7gYAQC6XQygUmrobhOOoQBFiIfz8/FBTU2PqbgAAZDIZRCKRqbtBOI4KFCEWQiQSobq62tTdAEAFiugHfQdFiIXw8vJCfX09mpub+3Sx7s8//9ynfsjlctTX18PT07NP6yGEjqAIsRB8Ph9RUVHIyckxaT9ycnIQHR3N+oxCQp5GBYoQC7Js2TJkZWVBo9GYpH2NRoOsrCzEx8ebpH1iWahAEWJBxGIxnJycUFRUZJL2i4qKIBQKIRaLTdI+sSxUoAixIDweD8nJycjMzIRSqTRq20qlEhkZGUhOTgaPxzNq28QyUYEixMKEh4djypQp2LNnj1HbzczMRFBQEKZOnWrUdonlogJFiAVau3YtiouLUVBQYJT2CgoKUFxcDBcXF1RWVhqlTWL56DQbQiyQo6MjUlNTsXr1agwYMAD+/v4Ga6u0tBTp6enYvHkzBg0ahGPHjuH27duIioqCtbW1wdollo/HMAxj6k4QQgyjsLAQ69evx8qVKxESEqL39efn5yM9PR0bNmxAcHAwAKCtrQ3Hjx9HS0sLFi1aRLc8IjqjAkWIhZNKpUhJSUFgYCASEhJgb2/f53UqlUrs3r0bJSUlSE1NfeasPYZhUFBQgLy8PMydOxfjxo3rc5uk/6ECRUg/oFAosHHjRuTl5eGtt95CUFAQrKzYfwWt0WhQVFSEjIwMhIWFYc2aNS98FEdtbS2+//57iEQizJgxg6b8CCtUoAjpR3Jzc7FlyxbI5XLExsYiMjISAoGgx+XkcjlycnKQlZUFgUCA5ORkhIeH96rN1tZWHD9+HO3t7Vi4cGGv2iMEoAJFCOd17cK9vfaIYRhIpVJIJBKcPXsWbm5u8PT0hKenJwQCAfh8PtRqNeRyOWQyGWQyGerr6xEdHY34+HiIxWLW1zkxDINffvkFhYWFWL58OZydnV/49yqVCpWVlZBKpeDxePDx8UFgYCCrNgn3UYEihKM0Gg3u378Pd3d37f/ZTtup1WrIZDKUl5ejvLwcTU1N6OjogK2tLYRCIUQiEUQiETw9PfVyb7179+7BxcUFNjY2z/2b9vZ2ZGRk4MqVK/Dw8IBQKIRMJkN4eDiWLVvW5z4Q7qACRQhHffPNNygrK8P8+fMRGRmpLU4Mw3D2Tg4Mw2DTpk1oampCZGQkoqOj0d7ejsrKSuzduxdJSUnw9fU1dTeJkdCFuoRwkEKhQFVVFaysrHDu3Dl8/fXX2mdB8Xg8cPVzp0Qiwb179zBv3jxER0cDAOzs7DB27FgAj54zRfoPKlCEcJCDgwPEYjHs7Ozg4+ODmzdvIjMzEydPnkRbWxt4PB5qampQX19v6q722rVr1/Dzzz9j9uzZ2mu21Go1AMDW1hZtbW0vnBoklocKFCEcZGVlhZiYGLi4uOCll17Cv//7v8PZ2RnZ2dnYuXMnLl26hC1btuDhw4cme/QGG0qlErt27UJwcDBCQkJgbW0NjUaj/d7r3LlzsLW1hZubm4l7SoyJChQhHDVkyBAIBAIcOnQIPj4+ePvttxEeHo67d+9i//79sLKygrOzs07XOxnblStX0N7ejpCQEAgEAjAMo+33rVu3UFhYCGdnZ/j4+Ji4p8SYzH/kEkKe0fUd08yZM1FVVQWZTAaBQIDFixdj/PjxUKlU6OzsxN69e3HlyhUT97ZnN27cgLOzs/aegV1HfY2Njbhw4QJkMhkWL14Me3t7zn6/RtijAkUIB3Wdpefj4wM/Pz9cunQJwKM39NzcXMyePRsxMTGoq6vDL7/8Ysqu9krXdVEtLS1gGAbW1tZoa2tDTk4Ozp8/j3nz5mH8+PGcPkORsEd3MyeEo7rerP39/ZGZmYmoqCgcOXIEQ4cOxZw5czBgwAD4+PjA1dXV1F3t0ejRo3H+/Hlcu3YNY8eOhb29Pfbu3YvKykpMnjwZCxcufO6yjxctlUoFPp8PlUoFW1tbY3WfGAhdB0WIBdi+fTvu3r2L+vp6vPPOO/Dx8eHckcbJkydx8uRJODs74+HDh3BwcEBERATmz58PPp//zIXIXYWJYRg8ePAAx48fx+3bt8Hn8zFgwADMmDEDQUFBJnxFpK+oQBHCYV1v0pcvX8b27dsxa9YsLFy4kLOnY9+6dQu1tbXg8XiYMGECHBwctGf0dVec1Go1CgoKcPToUdjb2yMoKAiOjo548OABzp8/jxkzZmDp0qUmfEWkL6hAEWIhlEol1Go1nJycTN0VvXm6MD39//Pnz+Onn37CgAED8N5778HFxQXAo9sllZWV4fjx45g9ezbCwsKM3nfSd3SSBCFm7OnPj919nmQYBhqNBvb29nBycrKos9ysrKzAMAzOnDkDtVr9RHG6fv06Lly4gPr6enR0dODYsWPYu3cvgEd3nxg/fjymTZtGD0zkMCpQhJix5uZmtLa2oqGhAcC/bmP0eBHqumZILpdr/8aSnDp1CgcPHsTPP/+sfd1tbW24fPkybt68iYCAAISHh2P06NGQSqXYuHEj2traYG9vj5kzZ0IkEpn4FRBd0Vl8hJihlpYWHD58GFKpFM3NzRg+fDiGDh2KuXPnwsPDA8C/prusrKxw//59fP3114iIiLC4BwNOmTIFbW1tEIlE2uLb0dGBy5cvY/LkyViyZIn2TMWRI0di+/btKCsrw+TJky0qh/6IjqAIMUPp6emora1FZGQk3nrrLQwePBgymQwff/wxjhw5AuBf01/AoyOt2tpa/P777xb3puzi4oJFixZpCzMAXL16Fe3t7ZgxY4b2eycAcHV1RWdnJ27evAmg+ylRwh10BEWImSkvL8ft27fxpz/9CSNGjAAATJ48Gb/++isKCgqQlZWFyspKvPnmmxg2bBgAwMvLC5s2bcKAAQNM2XWD6zp7z9raGjweD0Kh8InvpVpbW9HW1oaBAwcCsLzpzv6GjqAIMTNdF5l2nSquVqvB4/Hg6+uLxYsX49VXX4VSqYREItHerZxhGDg4OHDivnt90VVwhgwZAqFQiLq6Ou1RUtej5R0cHOj6JwtBR1CEmJmBAweioaEBtbW1cHd3116kyuPx4OjoiMjISNjZ2WHPnj3Iy8vD/Pnz+92RgpeXFzw9PbFnzx7cu3cPbW1tuHnzJiorK7Fs2bInpv0Id9F1UISYoR07duDXX3/FK6+8gpCQEG0Bevw6oP379+PXX3/Ff//3f/e7AtXlxIkTKC4uRktLCwYOHIi5c+fC399fL4+nJ6ZHBYoQM9TQ0ICtW7eitbUVkZGRCAoK0h4VdBWpM2fO4OLFi/jggw8gEAhM3GPTkcvlsLOz014L9jS6wSx3WfaENSEc5eLigpSUFPj5+eHgwYM4ePAgLl++DJVKBSsrKzQ0NOD69esQCoX9ujgBgEAggJ2d3XNPEFEqlaisrDRyr4g+0BEUIQamVqtRXV2N8vJyVFRUoLGxESqVCjY2NnB2doafnx9EIhG8vLy6nZq6desWvvrqKzAMAzs7O3h4eODevXtobW1FSkoKJ+5WbkoPHjzAoUOHMHr0aMTExPR5+q+v25P0HhUoQgyAYRiUlZVBIpEgOzsbbm5u8PLywpgxYyAUCsHn86FWq9HU1ISamhpUV1ejvr4e0dHRiI+Ph1gsfmJaSqPRoLCwEBUVFbh79y78/f3h6+sLT09PE75K7mhvb8fJkyfR2NiIRYsWsT6JQt/bk/QOFShC9Cw3NxdpaWlQKBSIjY1FZGRkr27g2tzcjOzsbGRlZUEgECA5ORnh4eFG6HH/wDAMioqKtA909PX17dVytD1NhwoUIXrS3NyMTZs2IT8/H4mJiQgKCtLpuiSNRoOioiJkZGRgypQpWLduHRwdHQ3Q4/7pt99+w7FjxzBu3DhERUU9984bhtieoaGhWLt2LW3PXqICRYgeSKVSpKSkIDAwEAkJCd2eTcaWUqnEnj17UFxcjNTUVIjFYj30lACPbjZ74sQJKBQKLFq06Jk7ntP2NA9UoAjpo8LCQqxfvx4rV65ESEiI3tdfUFCAbdu2YcOGDQgODtb7+vsrhmFQUFCAvLw8zJkzB97e3gBoe5oTKlCE9IFUKsXq1avx/vvvw9/f32DtlJaWIi0tDZs3b6ZP3npWV1eH77//Hn5+fnBzc8Of/vQn2p5mggoUITpSKBSIj49HYmKiQT5pP62goACZmZmQSCT0HYaeKZVKHDp0CN98843BjpyeRtuzZ3ShLiE62rhxIwIDA43yZgYAISEhCAgIwKZNm4zSXn9ib2+PqqoqhISE0PY0I1SgCNFBbm4u8vPzkZCQYNR2ly9fjry8POTm5hq1XUvXtT3feOMNo7ZL2/PFqEARwhLDMEhLS0NiYqJezu5iw97eHomJidiyZQs9jE9PaHuaLypQhLAklUqhUChM9syhoKAgyOVySKVSk7RvaWh7mi8qUISwdODAAcTGxprs4YBWVlaIjY2FRCIxSfuWhran+aICRQgLarUa2dnZiIyMNGk/IiMjcfbsWajVapP2g+toe5o3KlCEsFBdXQ03N7de3Yvtca2trZg1axaWLFkClUql/fn58+cxbtw47Nq1i9X6BAIB3NzcIJPJWC1HnkTb07xRgSKEhfLycnh5ebFezsHBAV988QUqKirw+eefAwDq6+uxbt06xMTE6HT2mKenJ8rLy1kvR/6Ftqd5owJFCAsVFRUYM2aMTsuKxWKkpKQgPT0dFy9exLp162BtbY1PPvlEp/XRG1rf0fY0b/Q0LUJYaGxsxNChQ3VePikpCRcuXMCKFSugUqmQmZmJQYMG6bQugUCAu3fv6twXQtvT3NERFCEsqFSqPj0llcfjYfHixejo6IBIJEJERITO6+Lz+ejo6NB5eULb09xRgSKEBRsbmz6daXX//n18/PHHEIvFKC8vxzfffKPzutRqNWxtbXVentD2NHdUoAhhwdnZGU1NTTotyzAMUlJSYGtri127diEpKQkbNmxARUWFTuuTy+XPPMeIsEPb07xRgSKEBT8/P9TU1Oi07I4dO3Dx4kVs2rQJQqEQKSkp8Pb2xgcffIC2tjbW65PJZBCJRDr1hTxC29O8UYEihAWRSITq6mrWy5WVlWHjxo344x//iNDQUACAra0tvvjiC9TW1uJvf/sb63XSG1rf0fY0b3QWHyEseHl5ob6+Hs3Nzawu7pwwYUK3pxB7eXmhrKyMdT/kcjnq6+vh6enJelnyL7Q9zRsdQRHCAp/PR1RUFHJyckzaj5ycHERHR/fpDDRC29PcUYEihKVly5YhKysLGo3GJO1rNBpkZWUhPj7eJO1bGtqe5osKFCEsicViODk5oaioyCTtFxUVQSgUQiwWm6R9S0Pb03xRgSKEJR6Ph+TkZGRmZkKpVBq1baVSiYyMDCQnJ4PH4xm1bUtF29N8UYEiRAfh4eGYMmUK9uzZY9R2d+/ejbCwMEydOtWo7Vo62p7miQoUITpau3YtiouLUVBQYJT2CgoKUFJSgjVr1hilvf6Gtqf5oQJFiI4cHR3xP//zP/jyyy9RWlpq0LZKS0uRnp6O1NRUODo6GrSt/srR0RGpqalIT083yvbctm0bbc8eUIEiREctLS0oKSlBUlIS/v73vxvsk3d+fj7+/ve/49NPP6Uv0g1MLBbj008/Nfj2/OKLLxAeHo4hQ4YYpA1LwWMYhjF1JwjhmtbWVuzbtw8+Pj6YNm0arl27hpSUFAQGBiIhIQH29vZ9bkOpVGL37t0oKSlBamoqFScjkkqlBt+ewKPrn5YtWwY3N7c+r98S0REUISy1trZi//798Pb2xrRp08Dj8SAWiyGRSGBnZ4eUlBRcvnxZ5+tqNBoNLl++jJSUFNjb20MikVBxMjJjbE+xWIyZM2fiwIEDePDggZ5fgWWgIyhCWFAqldi3bx/Gjh2LGTNmdHtqcG5uLrZs2QK5XI7Y2FhERkZCIBD0uG65XI6cnBxkZWVBIBAgOTkZ4eHhhngZhIXHt2dMTAyioqL0uj1LS0vx888/47XXXoOrq6shXgJnUYEipJeUSiX2798PT09PzJw584XXrTAMA6lUColEgrNnz8LNzQ2enp7w9PSEQCAAn8+HWq2GXC6HTCaDTCbDvXv3MGXKFKxYsQJisZiuizEjDMPgwoUL2Lp1K27evNmr7VlfX4/o6GjEx8f3uD27itTrr7+u8xN5LREVKEJ6QalU4sCBAxg9ejQiIyNZFQ+1Wg2ZTIby8nKUl5ejqakJHR0dsLW1hVAohEgkgkgkgkKhQE1NDZYsWWLAV0J0df78eahUKsyYMaNX29PT05PVvfWuXr2KCxcu4LXXXqMi9f9RgSKkB21tbdi/fz9GjRqFqKgogx3ZKJVKbN26Fe+99x49WdXMMAyD7du3Y/78+Rg2bJjB2ikpKcHFixfx+uuvw8XFxWDtcAWdJEHIC7S1teHAgQMYOXKkQYsTANjb22PYsGE6PZ+IGFZ9fT1UKhU8PDwM2k5AQADCw8Oxb98+NDQ0GLQtLqACRchztLW1QSKRYPjw4YiOjjbKd0J+fn46PzKcGE5lZSV8fX2NMgYCAwMRHh6O/fv3o7Gx0eDtmTMqUIR0o729HRKJBMOGDUNMTIzRTljw9vaGTCaDSqUySnukdyoqKuDn52e09gIDAxEWFoZ9+/b16yJFBYqQp7S3t+PAgQPw8PAwanECAAcHB3h4eNA0nxmpr69He3u7Qb976s6kSZMQGhqK/fv3o6mpyahtmwsqUIQ8puvIaejQoYiNjTXJqd40zWdejDm997SgoCCEhIRg3759/bJIUYEi5P9rb2/HwYMHMWTIEMyaNctk1yH5+PjQNJ8ZMfb03tOCg4MxefJk7Nu3D3K53GT9MAUqUITgX8XJzc0Ns2fPNulFsg4ODhg6dChkMpnJ+kAeefDgAdra2jB8+HCT9mPy5MkIDg7ud0WKChTp9zo6OvDdd9/B1dUVL730klncwcHX15em+cyAKaf3nhYSEoJJkyb1qyJFBYr0a13FycXFBXPmzDGLNyLg0TRfdXU1TfOZWEVFBXx9fU3dDa0pU6YgMDAQ+/fvR3Nzs6m7Y3BUoEi/pVKpcOjQIQiFQsydO9dsihMADBw4EO7u7qipqTF1V/qthw8forW11eTTe08LDQ3FxIkTsW/fPosvUlSgSL/UVZwEAoHZFacuNM1nWl3Te1ZW5vc2GRYWBn9/f+zfvx8KhcLU3TEY80ueEAPrKk6Ojo6YO3euWb4BAY+m+aqqqqBWq03dlX7J3Kb3njZ16lSIxWLs27fPYouUee6ZhBiISqXC4cOHMXDgQMybN89sixMAODo6YsiQITTNZwINDQ1oaWnBiBEjTN2VFwoPD8f48eOxf/9+tLS0mLo7eme+eycheqZWq3H48GHY29tj/vz5Zl2cutA0n2lUVlbCx8eHE2MkIiICfn5+2Ldvn8UVKfNPnxA96CpOAwYMwMsvv8yJNx7g0TTfjRs3aJrPyMx9eu9p06ZNg6+vL/bv34/W1lZTd0dvuLGXEtIHarUaR44cga2tLRYsWMCZ4gQATk5OGDx4MG7evGnqrvQbjY2NaG5uxsiRI03dFVamTZsGb29viypS3NlTCdGBWq3G0aNHYWNjw7ni1IWm+YyLS9N7j+PxeJg+fTrGjRtnMUWKW1uAEBa6ipO1tTUWLFgAa2trU3dJJ13TfJ2dnabuSr/Atem9x3UVqbFjx+LAgQNQKpWm7lKfUIEiFqmzsxPHjh2DtbU1Fi5cyNniBAACgQCurq50Np8RNDU1oampCaNGjTJ1V3TG4/EwY8YMeHp6Yv/+/ZwuUlSgiMXpKk48Ho/zxamLr68vKisrTd0Ni1dRUcHJ6b2n8Xg8zJw5E2PGjOH0kRS3twIhT+ns7MT3338PhmGwaNEiiyhOwKMCdf36dZrmM7Cuu0dYAh6Ph8jISIwaNQoSiQRtbW2m7hJrfFN3gPRfarUa1dXVKC8vR0VFBRobG6FSqWBjYwNnZ2f4+flBJBLBy8sLfH7PQ7WrOHV2dmLx4sUWU5yAR9N8gwYNQlVVFQDoLbP+ojdjbeTIkXj48CGnp/eexuPxEBUVhbNnz+LAgQNYtmwZBgwY0Ktl9b1/6tR/hmEYg6yZkG4wDIOysjJIJBJkZ2fDzc0NXl5eGDNmDIRCIfh8PtRqNZqamlBTU4Pq6mrU19cjOjoa8fHxEIvF3d43r7OzE//85z+hUqmwZMkSi3pz7sps69atuHLlCoYMGaKXzCwd27F248YN1NfXIzY21uJyYxgGZ86cwZ07d7Bs2TLY2dk99+8MsX/qigoUMZrc3FykpaVBoVAgNjYWkZGRcHJy6nG55uZmZGdnIysrCwKBAMnJyQgPD9f+XqPR4J///Cc6OjosrjgZKjNLR7k9i2EYZGVl4bfffuu2SJljZlSgiME1Nzdj06ZNyM/PR2JiIoKCgnT6Elqj0aCoqAgZGRkIDQ3F2rVr4eDggB9++AFtbW2Ii4uzmOJkyMwcHR0N0GPzQLm9GMMwOH36NH7//XfEx8fDzs7OrDOjAkUMSiqVIiUlBYGBgUhISIC9vX2f16lUKrFnzx4UFxdjwYIFGDRoEF555RWLKU6Gziw1NRVisVgPPTUvlFvvPF6k/P398R//8R9mmxkVKGIwhYWFWL9+PVauXImQkBC9r7+goABffvklNmzYgNDQUL2v3xSMkdm2bduwYcMGBAcH6339pkK5scMwDL7++mvs27cP7777rtlmRgWKGIRUKsXq1avx/vvvw9/f32DtlJaWIi0tDZs3b+b8p1vKTDeUG3tcyYwKFNE7hUKB+Ph4JCYmGuST2dMKCgqQmZkJiUTC2e8JKDPdUG7scSkzulCX6N3GjRsRGBholMEPACEhIQgICMCmTZuM0p4hUGa6odzY41JmVKCIXuXm5iI/Px8JCQlGbXf58uXIy8tDbm6uUdvVB8pMN5Qbe1zLjAoU0RuGYZCWlobExES9nA3Ehr29PRITE7FlyxZwadaaMtMN5cYeFzOjAkX0RiqVQqFQICgo6JnfrVq1CpMmTUJ7e3u3yyoUCkyYMAEpKSkAgDt37uC9995DQEAAAgIC8M477+DOnTsvbD8oKAhyuRxSqbTvL8ZIXpSZMXAxM4By04W+9s+TJ0/i3XffxfTp0zF+/HjExsYiNTUVCoXihe3rkhkVKKI3Bw4cQGxsbLcX+cXFxUEulyM7O7vbZX/88UcolUrExcVBqVRi+fLlqKqqQmpqKj777DPU1NQgISHhhQ9hs7KyQmxsLCQSid5ek6G9KDNj4GJmAOWmC33tn9u3b4e1tTXWrVuHb775BgkJCdizZw8SExOh0Wie274umVGBInqhVquRnZ2NyMjIbn8fFRUFFxcXHDlypNvfHz58GMOGDUNYWBj279+P27dvY+vWrZg9ezZmzZqFbdu2oa6uDvv27XthPyIjI3H27Fmo1eq+viSD6ymz3poxYwY2b96s8/JcygzQX259xaXc9Ll/btu2DWlpaVi0aBFCQ0ORlJSEDz/8EMXFxfjll19e2A+2mVGBInpRXV0NNze35967y9bWFi+//DLOnTuHhoaGJ353584d5OfnY/HixeDxeDhz5gwCAwMxZswY7d+MHDkSwcHByMrKemE/BAIB3NzcIJPJ+vyaDK2nzIyFS5kBuuXW2tqKWbNmYcmSJVCpVNqfnz9/HuPGjcOuXbtY94NLuelz/3R1dX1m+YkTJwIAfv/99xf2g21mVKCIXpSXl8PLy+uFfxMXFweVSoUffvjhiZ8fPXoUDMMgLi4OAHD9+nX4+Pg8s7y3tzdu3LjRY188PT1RXl7Oovem0ZvMjIUrmQG65ebg4IAvvvgCFRUV+PzzzwEA9fX1WLduHWJiYvDGG2/o1Beu5KbP/bM7+fn5AICxY8f22Bc2mVGBInpRUVHxxBFPdyZOnAhvb+9nphGOHj2KSZMmwdPTE8Cjx24LhcJnlhcKhWhqauqxL1x50+hNZk9jGAZqtfqJf8CjG3U+/jO2DzbkSmaAbrkBgFgsRkpKCtLT03Hx4kWsW7cO1tbW+OSTT3TuC1dy0+f++bS7d+/i888/R0REhPZI6kWoQBGja2xs7LaoPG3JkiUoKSnRHuKXlJSgqqrqhZ/O2BIIBL0qZKbW28wel5eXB19f3yf+1dXVYcuWLU/8jO0RAVcyA3TLrUtSUhKmT5+OFStW4MKFC0hNTcWgQYN07gtXcjPU/tnS0oJVq1aBz+fj008/7VVf2GRGBYrohUql6tXdxBcvXgwrKyscPnwYAHDkyBHY2tpi/vz52r953gB+3pHV0/h8Pjo6Olj03jR6m9njJkyYgCNHjjzxb8iQIVi2bNkTP/vf//1fVuvlSmaAbrl14fF4WLx4MTo6OiASiRAREdGnvnAlN33un13a2trwhz/8Abdv30ZGRgY8PDx61Rc2mVGBInphY2PTqzNz3N3dERERgWPHjqGjowPHjx9HTEzME4XH29sb169ff2bZGzduYNy4cT22oVarYWtry+4FmEBvM3uco6MjJk6c+MQ/GxsbuLu7P/Eztt/RcCUzQLfcuty/fx8ff/wxxGIxysvL8c033/SpL1zJTZ/7J/Co4L333nsoKyvDjh074Ovr2+u+sMmMChTRC2dn514ftsfFxaGurg6fffYZHj58+Mz0QWxsLIqLi3Hr1i3tz2pra1FYWIiYmJge1y+Xy3WeAjImNpkZGlcyA3TPjWEYpKSkwNbWFrt27UJSUhI2bNiAiooKnfvCldz0uX9qNBqsWbMGly5dwldffYVJkyax6gubzKhAEb3w8/NDTU1Nr/529uzZcHR0xM6dO+Hq6ooZM2Y88ftly5Zh+PDhWLVqFU6fPo2srCysWrUKHh4eeP3113tcv0wmg0gk0uVlGBWbzAyNK5kBuue2Y8cOXLx4EZs2bYJQKERKSgq8vb3xwQcfoK2tTae+cCU3fe6ff/nLX3DixAm8/fbbcHBwwJUrV7T/fvvttx7XzyYzKlBEL0QiEaqrq3v1twMGDMC8efPAMAwWLlz4zNy4g4MD9uzZA09PT6xbtw5r1qzBiBEjsHv3bgwcOLDH9XPlTcPX1xdVVVWm7gYA7mQGsBtrXcrKyrBx40b88Y9/1D7c0tbWFl988QVqa2vxt7/9Tae+cCU3fe6f586dAwB8+eWXWLp06RP/enOXCDaZWcYzsonJeXl5ob6+Hs3Nzb26gPKTTz554em9w4YNw5dffsm6H3K5HPX19c89JdbUNBoN6urqUFFRgfLycty/f7/XmT3Pzz//3Kc+mXtmT2M71oBHJ5d0d2qzl5cXysrKdOoHl3LT5/7Zl/HGNjM6giJ6wefzERUVhZycHJP2IycnB9HR0Tqf5WUIDMOgtrYWWVlZ+Oqrr5CVlQVHR0csX74cMTExlBlLNNbY42pmVKCI3ixbtgxZWVkvvGGkIWk0GmRlZSE+Pt4k7T/u8aL05Zdf4tSpU3BwcMBrr72GpKQkTJ06FYMGDaLMdES5scfFzMy/9BPOEIvFcHJyQlFRESZPnmz09ouKiiAUCiEWi43eNvCoKN25cwcVFRWorKyEra0tRCIRli1bBjc3t26X6e+Z6YpyY4+LmdERFNEbHo+H5ORkZGZmQqlUGrVtpVKJjIwMJCcng8fjGa3drqJ09uxZbN26FSdOnICtrS1effVVrFixAhEREc8tTkD/zEwfKDf2uJgZj+HSIyEJJ/z1r39Fe3s7VqxYYbQ209PTYW9vjw8//NDgbTEMg7t372qPlKytreHr6wuRSAQ3Nzed3rQsPTNDodzY41JmVKCI3ikUCsTHxyMxMREhISEGb6+goACZmZmQSCRwdHQ0SBsMw+D3339HRUUFKioqYGVlBT8/P/j5+WHw4MF9/iRtiZkZgyly27lzJw4dOsTZ3Lg01ug7KKJ3jo6OSE1NxerVqzFgwAD4+/sbrK3S0lKkp6dj8+bNen/DYBgG9+7d0xYl4NEFj0uWLMGQIUP0Or1jKZkZm7Fz27ZtG2bMmAGpVIopU6ZwaoqvC5fGGh1BEYMpLCzE+vXrsXLlSoN8UsvPz0d6ejo2bNiA4OBgvayTYRjcv39fW5Q0Go32SMnd3d3gb0hczMwcGDM3b29vHDp0CB4eHpg9e7bJHjvfV1wYa1SgiEFJpVKkpKQgMDAQCQkJsLe37/M6lUoldu/ejZKSEqSmpurlTKrHi5JardYWpaFDhxr9UzJXMjM3xsytvb0d33//PRiGwaJFi2BnZ9fntkzB3McaFShicAqFAhs3bkReXh7eeustBAUF6fSpU6PRoKioCBkZGQgLC8OaNWv6NEVVX1+vLUodHR3aouTh4WHyqRtzzczcGTM3jUaD06dPo66uDkuXLoVAINDXyzAqcx5rVKCI0eTm5mLLli2Qy+WIjY1FZGRkr3ZquVyOnJwcZGVlQSAQIDk5GeHh4Tr14cGDB9qi1NbWpi1Kw4YNM3lRAh5NMT7eD3PIjIuMlRvDMMjPz0dhYSFeeeUVuLu76/NlGJU5jjUqUMSoGIaBVCqFRCLB2bNn4ebmBk9PT3h6ekIgEIDP50OtVkMul0Mmk0Emk6G+vh7R0dGIj4+HWCxmXUgePnyoLUpKpRK+vr7w8/PD8OHDzaIo9cQUmVmCp3MTCoUYN24cxo4dq/fcKisrcerUKcybNw9jx4418CszHHMba1SgiMmo1WrIZDKUl5ejvLwcTU1N6OjogK2tLYRCIUQiEUQiETw9PVnf76yhoUFblFpaWrRFacSIEWb7Zp2bm4vr169j0aJFcHZ2hkajeWaqxZCZWTK1Wo1PP/0UXl5euHXrlkFyq6urw9GjRzF16lQEBQXp+RUYnzmMNSpQxKQ6OzthbW2tl3U1NjZqi5JCoYCPj4+2KHHhTKvPPvsMN27cQFhYGF599dVePVqE9E5HRwfS0tLwpz/9yaBjobGxEQcPHsTYsWMRFRVlth+GuII+YhGjU6lUyMvLQ3FxMQBgzJgxCAkJ0Wn+vqmpSVuU5HI5fHx8EBUVhZEjR3KiKHVRKpWoqqrC1KlTceXKFVy/fh1Lly5FYGAgeDzeM99NEXYePHiAQYMGGXxMODs7Y/ny5Thy5AiOHj2Kl19+GTY2NgZt05LRERQxmvb2dpw+fRq5ubmwsrLC+PHj0d7ejry8PHh5eWH58uUYNmxYj+uRy+XaotTY2Kg9Uho1ahSnitLjzp07h9OnT+Ojjz7C1atX8d1330GlUmHp0qXaB+wR3ZWWlqKmpgYLFiwwSntqtRo//vgjHj58iFdeeYWOhnVER1DEaNRqNe7du4ewsDBMnz4dTk5O4PP5CAgIwLfffotr165h2LBh3R4tMAyDsrIyFBcXo6GhAd7e3pg+fTpGjRqltylCU8rNzcWECRPA5/MRFBQEDw8PHDp0CBkZGbh+/ToWLlzI2dOYzcGDBw/g6upqtPb4fD7mz5+PixcvYteuXVi6dOkLbxpMukcFihjNgAEDMH/+fAwaNAg2Njba59IMGzYM7e3t2h24u6mszs5O2NvbIyIiAqNHj7aIotRFpVLh1q1beOWVV7Q/8/DwwOuvv47vv/8eFy9ehFwuR1xcHIYOHWrCnnJXfX09Jk6caNQ2eTwepk2bBqFQiH379mHhwoUYPXq0UfvAddycDyFmR61W9/g31tbWcHd3187Jd03HFRUVYfz48fD29n7usnw+H+PGjYOXl5dFFScAKC4uhpOTE4YPH/7Ez11dXZGUlIQ333wT1dXV+Oqrr3R+PHl/V19fb7IjGH9/fyxcuBDff/89SktLTdIHrqLvoAhrXVNwtbW1OHXqFG7duoURI0ZgwoQJCAsL69XyarUaFy5cwKlTp9DQ0ICRI0di3LhxiIyMhLu7e786KUAqlaKtre2Z+5V1ZaDRaHDlyhUcOnQILS0tmD9/PmbPnm2i3nKPsc7g60l9fT2+++47iMViTJs2rd+M776gIyjCGo/Hw4MHD7Bz504oFApMmzYNra2tyMjIwE8//YT29nYAj95gn7e8XC7Hb7/9hqlTp+Kdd95BcHAwLl++jK+++goqlapf7bxisbjb62a6MrCyskJwcDD+8z//E6NGjcKdO3eM3UVOe/jwIVxcXEx+Ao2bmxveeOMNyGQyHD9+vFezDv0dHUERnezcuRMymQzJycna08MlEgmKioqwYMECREREdHuh6eO6ft91pJCXl4eMjAwkJiYiLCysXx1FvQjDMGAYBlZWVrh//z74fD5cXFxM3S3OKCsrQ3V1NRYuXGjqrgB49J3jDz/8AKVSibi4OAwYMMDUXTJbdARFWGMYBr/99ht8fX3h7u6uPdkhOjoabm5uyM7OBoAeP7E+/fvRo0fD3t4ed+/eBdD9yRL9EY/H02Y1ePBgKk4sGfsMvp7Y2Nhg0aJFGDp0KHbt2oXGxkZTd8lsUYEirHV2dsLKyko7lddVSNzc3BAQEIDa2lrcvn0bwPOn+R7X9T1LRUUFWltb4ePjY7jOcxxNeLBnyhMknsfKygrR0dEIDg7G7t27UVdXZ+oumSUqUIQ1Pp+PgQMHQqlUoqWlRVtgAGhvxHnp0iUA/3pD7fr945RKJZqamlBfX4+LFy/i559/xpw5cyASiYz3YjhCpVIBoKNKXZhjgeoSFBSEuXPn4tChQ6isrDR1d8wOXQdFdBIREYGDBw/i3r178PT01P588ODBGDVqFK5fvw7gX9N4T0/nXb16FQUFBZDL5airqwOPx0NERASio6P75ZtwVyF//LV3fQfX0tICiUQCPz8/TJ061VRd5CSVSgWFQgFnZ2dTd+W5xo4di2XLluHQoUNobGzk7KPkDYGOoIhOxGIxGhsbIZPJAPyrADk5OUEgEECj0Twxty6VSrFjxw7cvHkTADB8+HB4eXlh1KhRSEpKQmpqKhYvXtxv75bA4/G0b0oajeaJE0RqamqQl5dH19Do4OHDh3B2djb7a+fc3d2xfPlySKVSnDp1qtsZh/6IjqBIj7o7m27AgAHw8fHB1atXMWnSJLi4uGjvTD5gwAAwDPPEY7B///13FBQUYPz48Rg9ejRcXV0RFRVl7JdiVtra2lBfX4/a2lpoNBp4enrCw8PjiaNNhmEgFovxf//3f/T9kw7MeXrvaQKBAAkJCTh27BgOHTqEhQsXcvZR8vpCBYp0q6OjAzdu3EBlZSXGjx8PHx+fZ4rUokWLsG3bNly4cAELFiyAtbU1mpubIZPJIBAIYG9vr/3b8PBwDBkyBGKx2NgvxSxdv34dJ0+exLVr12BnZwc7OzsoFAoIBAIEBQVh2rRp2qf8Mgxj1lNU5szczuDriZ2dHV555RWcPn0ae/fuxSuvvNJvZxUAug6KPKajowNVVVWoqKhATU0NRowYAT8/P/j6+sLW1vaZv9doNDh8+DCysrIwffp0TJgwAdeuXUNZWRn+8Ic/YMyYMcZ/ERzx0UcfwdnZGQsWLICdnR3u3buHBw8eQCaT4caNG2AYBmFhYZg3b94ThZ6wc/jwYYwfPx5+fn6m7gorlvQo+b6gI6h+TqVSaYuSTCbD8OHD4efnhzlz5vT4xmhlZYWlS5fCzs4OlZWVKC4uhlAoRFxcHN0U8wWqqqogl8uRnJysnX4aMWIEAKClpQW3bt1Cfn4+srOzoVQq8dprr9HTcXXEpSm+x/F4PISGhkIoFEIikXD+UfK6olHfD6lUKlRXV2uLkoeHB/z8/PDSSy/p9Gl9wYIFiImJAY/Ho0/7vfDbb7/B0dERnZ2dAB7daNfKygpWVlYYOHCg9lHaIpEI3377LaZNm0ZHozpQq9WQy+WcvrDZz88Pjo6OOHr0KCIiIjBp0iRTd8moqED1E2q1WluUqqurMXToUPj5+WHWrFlwcHDo8/r1sY7+QiQS4ciRIygtLYW7u/sTR0cajQYajQZ8Ph8ikQjDhw/HlStXqEDpgCtn8PVkxIgRSEhIwHfffYfGxkZERkb2m9PQqUBZMLVaDZlMhoqKClRVVcHd3R1+fn6IiYmhJ3ya0KBBgzB16lQcOnQIN2/eRHh4OEaPHg0HBwftkRTw6JY4HR0ddIKEjrg6vdcdFxcXLF++HIcPH+5Xj5LvtydJdB1RlJeXax8drlKpYGNjA2dnZ/j5+UEkEsHLy4tT8/9qtRo1NTWoqKjAjRs3MGTIEPj5+cHHxweOjo6m7l6/9Lyx1tjYiIaGBgwcOBDjx49HcHAwhg4diiFDhsDW1hbHjh2DVCrFRx99RB8onuNF+7FSqcSIESOwcOFCzu3Hz6NWq3Hy5Ek0NjYiLi5Op3HBpfe+flWguh4bLpFIkJ2dDTc3N3h5eWHMmDEQCoXg8/lQq9VoampCTU0NqqurUV9fj+joaMTHx0MsFpvloXVnZ+cTRcnNzU179h0VJdNgM9aqqqpw/fp1PHz4EEOHDsXgwYNhZWUFDw8PzJs3D1OmTDH1yzErbLKtrq5GTU0NJ/bj3mIYBhcuXMC1a9ewdOnSXp1Gz9X3vn5ToHJzc5GWlgaFQoHY2FhERkbCycmpx+Wam5uRnZ2NrKwsCAQCJCcnIzw83Ag9frHOzk7cvHkTFRUVuH79OlxdXbVFqTevixhOX8faqVOn4ODggBUrVmDOnDlG6DF3WNp+3BdXr17FuXPnsGjRIowaNeq5f8flzCy+QDU3N2PTpk3Iz89HYmIigoKCdHpwmUajQVFRETIyMhAaGoq1a9ca/eiks7MTt27d0hYlFxcXbVHqzxfzmQtLGmvmhrLtXk1NDf75z38iOjr6mYvgLSEziy5QUqkUKSkpCAwMREJCgl5OgVYqldizZw+Ki4uRmppq8DsjaDQabVH69ddftXPEfn5+VJTMiCWMNXNF2b5Y16PkJ0yYgIiICPB4PIvJzGILVGFhIdavX4+VK1ciJCRE7+svKCjAtm3bsGHDBgQHB+t13RqNBrdv39YWJaFQCF9fX/j5+UEoFOq1LdJ3XB5r5o6y7R2FQoHDhw/D1dUVQ4YMwZ///GeLyMwiC5RUKsXq1avx/vvvw9/f32DtlJaWIi0tDZs3b+7zpwmNRoPa2lpUVFSgsrISAoFAO31HpxmbLy6ONa6gbNlRqVT4+uuvcfDgQaxZs8YiMrO4AqVQKBAfH4/ExESDfHp4WkFBATIzMyGRSFjPy2o0GtTV1WmLkqOjo7Yocfnq9/6CS2ONayhb9iwxM4srUH/961/R3t6OFStWGK3N9PR02Nvb48MPP+zxbxmGeaIoOTg4aIvSoEGDjNBboi/mPta4jLJlzxIzs6gHFubm5iI/Px8JCQlGbXf58uXIy8tDbm5ut79nGAa1tbXIysrCl19+qT2N+LXXXkNSUhKmTp1KxYljzHWsWQLKlj1LzcxijqAYhsG//du/YcmSJZg8ebLR2798+TKOHj2KPXv2aJ/hc+fOHe2Rkp2dnfZIyVJuv9JfmdtYsySULXuWnJnFHEFJpVIoFAoEBQWZpP2goCDI5XKcO3cOZ8+exdatW3HixAnY2dnh1Vdfxdtvv42IiAgqThbAXMaaVCo1SfuG9KJsV61ahUmTJqG9vb3bZRUKBSZMmICUlBQAwJ07d/Dee+8hICAAAQEBeOedd3Dnzp0Xts/FbPWV2W+//YaPPvoIS5cuhVgsxtixY1FbW9tj+4bMzGIK1IEDBxAbG6vThWj6YGVlhZiYGGzfvh02NjZYunQpVqxYgWnTpmHw4MEm6RMxDHMYa7GxsZBIJCZp35BelG1cXBzkcjmys7O7XfbHH3+EUqlEXFwclEolli9fjqqqKqSmpuKzzz5DTU0NEhIS0Nra+tz2uZitvjK7efMmTpw4AYFAwOokC0NmZhEFSq1WIzs7G5GRkSbtR1RUFGpqajB16lQMHjyYM1MEpPfMZaxFRkbi7NmzUKvVJu2HPvWUbVRUFFxcXHDkyJFuf3/48GEMGzYMYWFh2L9/P27fvo2tW7di9uzZmDVrFrZt24a6ujrs27fvhf3gUrb6zGzKlCnIz8/Hzp07MXfuXFb9MFRmFlGgqqur4ebmxvoedK2trZg1axaWLFkClUql/fn58+cxbtw47Nq1i9X6BAIB3NzcIJPJWC1HuIPGmuH0lK2trS1efvllnDt3Dg0NDU/87s6dO8jPz8fixYvB4/Fw5swZBAYGPvEcrZEjRyI4OBhZWVkv7AeXstVnZn2ZETBUZhZRoMrLy+Hl5cV6OQcHB3zxxReoqKjA559/DuDRbUPWrVuHmJgYvPHGG6zX6enpifLyctbLEW6gsWY4vck2Li4OKpUKP/zwwxM/P3r0KBiGQVxcHADg+vXr8PHxeWZ5b29v3Lhxo8e+cCVbfWbWV4bIzCIKVEVFhc5PHBWLxUhJSUF6ejouXryIdevWwdraGp988olO6+PKwCa6obFmOL3JduLEifD29n5myuro0aOYNGkSPD09AQBNTU3d3hZMKBSiqampx75wJVt9ZtZXVKCeo7GxsU/3qEtKSsL06dOxYsUKXLhwAampqTpflyQQCHq1AxBuorFmOL3NdsmSJSgpKdFOJ5WUlKCqqkpvRwIAd7K19MwsokCpVKo+PfmRx+Nh8eLF6OjogEgkQkREhM7r4vP56Ojo0Hl5Yt5orBlOb7NdvHgxrKyscPjwYQDAkSNHYGtri/nz52v/5nlvls87snoaV7LVZ2Z9ZYjMLKJA2djY9Onskfv37+Pjjz+GWCxGeXk5vvnmG53XpVarYWtrq/PyxLzRWDOc3mbr7u6OiIgIHDt2DB0dHTh+/DhiYmKeKDze3t64fv36M8veuHED48aN67ENrmSrz8z6yhCZWUSBcnZ21vnQkmEYpKSkwNbWFrt27UJSUhI2bNiAiooKndYnl8vpkRgWjMaa4bDJNi4uDnV1dfjss8/w8OHDZ6aqYmNjUVxcjFu3bml/Vltbi8LCQsTExPS4fq5kq8/M+soQmVlEgfLz80NNTY1Oy+7YsQMXL17Epk2bIBQKkZKSAm9vb3zwwQdoa2tjvT6ZTAaRSKRTX4j5o7FmOGyynT17NhwdHbFz5064urpixowZT/x+2bJlGD58OFatWoXTp08jKysLq1atgoeHB15//fUe18+VbPWZGQCcPHkSJ0+eRFlZGQDg3LlzOHnyJPLy8npcvyEys4gCJRKJUF1dzXq5srIybNy4EX/84x8RGhoK4NF1A1988QVqa2vxt7/9jfU6uTKwiW5orBkOm2wHDBiAefPmgWEYLFy48JnvYRwcHLBnzx54enpi3bp1WLNmDUaMGIHdu3dj4MCBPa6fK9nqMzMASE5ORnJyMvbu3QsA+PDDD5GcnIzNmzf3uH5DZGYRN4tVq9WIjIzEP/7xD9YXUOqTXC5HcnIycnJy+vRFOjFfNNYMh7Jlz9Izs4gjKD6fj6ioKOTk5Ji0Hzk5OYiOjjb7QU10R2PNcChb9iw9M4soUMCjOeesrCxoNBqTtK/RaJCVlYX4+HiTtE+Mh8aa4VC27FlyZhZToMRiMZycnFBUVGSS9ouKiiAUCiEWi03SPjEeGmuGQ9myZ8mZWUyB4vF4SE5ORmZmJpRKpVHbViqVyMjIQHJyMt3BvB+gsWY4lC17lpyZxRQoAAgPD8eUKVOwZ88eo7a7e/duhIWFYerUqUZtl5gOjTXDoWzZs9TMLKpAAcDatWtRXFyMgoICo7RXUFCAkpISrFmzxijtEfNBY81wKFv2LDEziytQjo6OSE1NRXp6OkpLSw3aVmlpKdLT05GamgpHR0eDtkXMD401w6Fs2bPEzCziOqjuFBYWYv369Vi5ciWrxxf3Vn5+PtLT07FhwwYEBwfrff2EO2isGQ5ly54lZWaxBQoApFIpUlJS4O/vjzfffBP29vZ9XqdSqcS3336L0tJSpKamcupsH2I4UqkU69atw8SJE/U61nbt2oWrV6/267HWtR8HBgYiISFBb9nu3r0bJSUlFpmtpWRmcVN8jxOLxfjzn/+Me/fuISUlBZcvX9b5WgGNRoPLly9j3bp1uH37NrZs2WJxg5robvz48XjttdfQ2tqqt7G2du1a3LlzB3v37u3XY00sFkMikcDOzk5v2aakpMDe3h4SicQis7WUzCz6CKqhoQG7du3C8uXLUVFRgS1btkAulyM2NhaRkZEQCAQ9rkMulyMnJwdZWVkQCARITk6Go6MjCgsL8eabb8La2toIr4SYu+LiYpSUlOCNN97AL7/8opex9t577+HevXsYNGgQZs6caYRXYf5yc3P1th+Hh4cbocemx+XMLLZAMQyDvXv3wsfHRzsPyzAMpFIpJBIJzp49Czc3N3h6esLT0xMCgQB8Ph9qtRpyuRwymQwymQz19fWIjo5GfHw8xGIxeDweGIbBwYMHMWLEiH4zyMnzyeVyZGRk4PXXX8fgwYMB6G+stbS0YOfOnVi6dCk8PDxM/ErNg76y7U+4mpnFFqjCwkJcu3YNCQkJsLJ6diZTrVZDJpOhvLwc5eXlaGpqQkdHB2xtbSEUCiESiSASieDp6dnt/aW6e1Mi/U9vPqz0daxdu3YNly5dQmJiIifuD2dMfc22P+JSZhZZoBobG/Htt98iISEBrq6uBmvn8Wmd7oogsXxXr15FUVER3njjDYNN9zIMgyNHjsDNza3bZ/gQYqks7l2VYRicPHkSoaGhBi1OABAQEAA7Ozvk5+cbtB1inpqbm3Hu3DnMmzfPoN9F8ng8zJ49GyUlJbh7967B2iHE3FhcgSouLoZKpTLI+f9P4/F4mDt3LvLz81FfX2/w9oj5YBgGP/30EyZNmoQhQ4YYvD1HR0dERUXhxIkT6OzsNHh7hJgDiypQTU1NOH/+PObNm2e0KTehUIjp06fjxIkTJrvdPTE+qVQKuVxu1Pu2icViCAQCXLp0yWhtEmJKFlOguqb2QkJC4ObmZtS2AwMDYWNjY7R7YBHTUigUyM7ONvjU3tN4PB5eeuklXLlyBffu3TNau4SYisUUqJKSErS3tyM0NNTobXdN9eXl5eHBgwdGb58YD8MwOHXqFAICAjB06FCjt+/k5ISZM2fSVB/pFyyiQMnlcvz8889Gndp7mrOzMyIiInDy5Ema6rNg5eXlaGhoMOn1b/7+/nBwcEBeXp7J+kCIMXC+QDEMgx9//BHBwcEmvx4pKCgIPB4PhYWFJu0HMYyWlhacPXsW8+bNM+n1ITweD3PmzMHly5dx//59k/WDEEPjfIEqKytDS0sLwsLCTN0V7VTfpUuX0NDQYOruED07ffo0JkyYYBZ3dBAIBJgxYwadnEMsGqcLVHNzs0m+rH6RQYMGISwsDCdPnoQFXgPdb1VUVOD+/fuYNm2aqbuiRdfhEUvH2QL1+HUo7u7upu7OEyZPnozOzk4UFRWZuitED1pbW5GVlWXyqb2n0XV4xNJxtkBdu3YNcrncLG/WamVlhXnz5uHixYtobGw0dXdIH2VlZWH8+PEYPny4qbvyDKFQiGnTptHJOcQicbJAKRQK7ZfV5jK19zRXV1eEhobSVB/H/frrr7h79y6mT59u6q4816RJk8Dn83H58mVTd4UQveJcgTL1dShshISEQKVSobi42NRdITpQKpU4ffo05s2bBxsbG1N357m6zur75Zdf8PDhQ1N3hxC94VyBKi8vx8OHD81yau9pXVN958+fR1NTk6m7Q1g6c+YMfH19MWLECFN3pUcuLi4IDw+ns/qIReFUgWppacGZM2fM7svqF3Fzc0NISAh+/PFHmurjkBs3bqCuro5Tj7cIDg4Gj8ejk3OIxeBUgTp9+jT8/f0xbNgwU3eFldDQULS1teHq1aum7grphba2Nvz000+YO3cubG1tTd2dXus6qy83N5euwyMWgTMFyhyvQ+mtrqm+c+fOQS6Xm7o7pAdnzpyBj48PRo0aZequsEbX4RFLwokC1XUdyty5czkztfe0wYMHIzg4GD/99BO9cZixqqoq3L59GzNnzjR1V3TWdR3elStXTN0VQvqEEwWq6zoULnxZ/SJhYWFQKBQoKyszdVdIN9rb2zk5tfe0riP2Cxcu0HV4hNPMvkBx4TqU3rK2tsa8efOQk5OD5uZmU3eHPOXs2bMYO3YsRo8ebequ9FnXdXh0cg7hMrMuUF3XocydO9esr0Nhw93dHYGBgTTVZ2ZkMhlqamoQGRlp6q7oTUhICDo6OlBSUmLqrhCiE4N/oaNWq1FdXY3y8nJUVFSgsbERKpUKNjY2cHZ2hp+fH0QiEby8vJ75fqnrOpSRI0cauptGFR4ejszMTFy7dg1isfiZ3/cls/6qL5m1t7fjxx9/xJw5c2BnZ2eiV6B/VlZWmDt3Lvbt2wdPT08IhcJn/obGGjFnPMYAH+MZhkFZWRkkEgmys7Ph5uYGLy8vjBkzBkKhEHw+H2q1Gk1NTaipqUF1dTXq6+sRHR2N+Ph4iMViVFVV4cyZM0hKSuL09wHPc/fuXRw8eBBJSUlwdHTUS2Y8Hs/UL8uo9JXZTz/9BI1Gg7lz55r6JRnEpUuXcOvWLcTHx4PH49FYI5yh9wKVm5uLtLQ0KBQKxMbGIjIyEk5OTj0u1/XojKysLDg5OWHMmDF45513OHmqb2+dO3cODx48gLu7O7Zs2dKnzAQCAZKTkzlxhw190Mc4EwgEiI+Px507d/D2229b1NHT4zo7O7Fr1y5MmjQJLS0tesmtP401Yjp6K1DNzc3YtGkT8vPzkZiYiKCgIJ0ev67RaFBUVIQdO3YgPDwca9euhaOjoz66aHYaGhqwevVq3L9/H0lJSX3OLCMjA6GhoRadmb7H2fbt2xEQEIC//OUvFpsZAFRXV+O//uu/0NjYiLfeeovGGuEEvRQoqVSKlJQUBAYGIiEhAfb29n3umFKpxJ49e1BcXIzU1NRuv6vhsq7MAgICsHz5csqsF2ic6aYrt4kTJ+KNN96g3Ahn9LlAFRYWYv369Vi5ciVCQkL01S+tgoICbNu2DRs2bEBwcLDe128KlBl7lJluKDfCZX0qUFKpFKtXr8b7778Pf39/ffbrCaWlpUhLS8PmzZs5/0mNMmOPMtMN5Ua4TucCpVAoEB8fj8TERIN8MntaQUEBMjMzIZFIODvnTZmxR5nphnIjlkDnC3U3btyIwMBAowx+4NFFhwEBAdi0aZNR2jMEyow9ykw3lBuxBDoVqNzcXOTn5yMhIUHf/Xmh5cuXIy8vD7m5uUZtVx8oM/YoM91QbsRSsC5QDMMgLS0NiYmJejkbiA17e3skJiZiy5YtnLpNEGXGHmWmG8qNWBLWBUoqlUKhUCAoKOiJn69atQqTJk1Ce3t7t8spFApMmDABKSkp+O233/DRRx9h6dKlEIvFGDt2LGpra3vVflBQEORyOaRSKduum8zzMgPY5Xby5Em8++67mD59OsaPH4/Y2FikpqZCoVC8sP3+nJmuY42LmQH62T91HWcAd3Mj5ol1gTpw4ABiY2OfucgvLi4Ocrkc2dnZ3S73448/QqlUIi4uDjdv3sSJEycgEAhYz5FbWVkhNjYWEomEbddN5nmZAexy2759O6ytrbFu3Tp88803SEhIwJ49e5CYmAiNRvPc9vtzZrqONS5mBuhn/9R1nAHczY2YJ1YFSq1WIzs7u9s7PkdFRcHFxQVHjhzpdtnDhw9j2LBhCAsLw5QpU5Cfn4+dO3fqdP+zyMhInD17Fmq1mvWyxvaizAB2uW3btg1paWlYtGgRQkNDkZSUhA8//BDFxcX45ZdfXtiP/ppZX8YalzID9Ld/9mWcAdzLjZgvVgWquroabm5u3d67y9bWFi+//DLOnTuHhoaGJ353584d5OfnY/HixeDxeDrdYuVxAoEAbm5ukMlkfVqPMbwoM4Bdbq6urs8sP3HiRADA77///sJ+9NfM+jLWuJQZoL/9sy/jDOBebsR8sdp7y8vL4eXl9dzfx8XFQaVS4Ycffnji50ePHgXDMIiLi9Otl93w9PREeXm53tZnKD1lBvQtt/z8fADA2LFje+wLZcYeVzIDDLt/shlnALdyI+aLVYGqqKjAmDFjnvv7iRMnwtvb+5lphKNHj2LSpEnw9PTUqZPd4coO0FNmgO653b17F59//jkiIiK0n3BfhDJjjyuZAYbbP9mOM4BbuRHzxapANTY2dvvQs8ctWbIEJSUl2sP7kpISVFVV6fXoCXg0jdDU1KTXdRpCbzID2OfW0tKCVatWgc/n49NPP+1VX/p7ZrrgSmaAYfZPXcYZwK3ciPliVaBUKlWPT9VcvHgxrKyscPjwYQDAkSNHYGtri/nz5+vey27w+Xx0dHTodZ2G0JvMAHa5tbW14Q9/+ANu376NjIwMeHh49Kov/TkzXXElM0D/+6eu4wzgVm7EfLEqUDY2Nj2emePu7o6IiAgcO3YMHR0dOH78OGJiYnr1iZgNtVrNiSft9iYzoPe5qVQqvPfeeygrK8OOHTvg6+vb677018z6giuZAfrdP/syzgBu5UbMF6sC5ezs3KvD9ri4ONTV1eGzzz7Dw4cP9T69BwByuVzvRc8QepsZ0HNuGo0Ga9aswaVLl/DVV19h0qRJrPrSHzPrK65kBuhv/+zrOAO4lRsxXz3PozzGz88P586d6/HvZs+eDUdHR+zcuROurq6YMWPGM39z8uRJAEBZWRmAR48/HzRoEAYNGoTQ0NAe25DJZIiKimLTfZPobWZAz7n95S9/wYkTJ/Duu+/CwcEBV65c0f5u6NChPU7B9MfMgL6NNa5kBuhv/+zrOAO4lRsxX6yOoEQiEaqrq3v8uwEDBmDevHlgGAYLFy7sdl48OTkZycnJ2Lt3LwDgww8/RHJyMjZv3tyrvshkMohEIjbdN4neZgb0nFvXm8+XX36JpUuXPvGvN1fu98fMgL6NNa5kBuhv/+zrOAO4lRsxX6yeB6VWqxEZGYl//OMfz72I0hjkcjmSk5ORk5PTqy/TTYkyY48y0w3lRiwNqyMoPp+PqKgo5OTkGKg7vZOTk4Po6GhODH7KjD3KTDeUG7E0rO8Ds2zZMmRlZfV400hD0Wg0yMrKQnx8vEna1wVlxh5lphvKjVgS1gVKLBbDyckJRUVFhuhPj4qKiiAUCiEWi03Svi4oM/YoM91QbsSSsC5QPB4PycnJyMzMhFKpNESfnkupVCIjIwPJycng8XhGbbsvKDP2KDPdUG7Ekuh0q+fw8HBMmTIFe/bs0Xd/Xmj37t0ICwvD1KlTjdquPlBm7FFmuqHciKXQ+VkEa9euRXFxMQoKCvTZn+cqKChASUkJ1qxZY5T2DIEyY48y0w3lRiyBzgXK0dERqampSE9PR2lpqT779IzS0lKkp6cjNTUVjo6OBm3LkCgz9igz3VBuxBKwug6qO4WFhVi/fj1WrlzJ+vHtvZGfn4/09HRs2LABwcHBel+/KVBm7FFmuqHcCJf1uUABgFQqRUpKCgIDA5GQkAB7e/s+d0ypVGL37t0oKSlBamqqxZ0VRJmxR5nphnIjXNW3Z6//f2KxGBKJBHZ2dkhJScHly5d1vg5Do9Hg8uXLSElJgb29PSQSiUUOfsqMPcpMN5Qb4Sq9HEE9Ljc3F1u2bIFcLkdsbCwiIyMhEAh6XE4ulyMnJwdZWVkQCARITk5GeHi4Prtmtigz9igz3VBuhEv0XqAIIYQQfdDLFB8hhBCib1SgCCGEmCUqUIQQQswSFShCCCFmiQoUIYQQs/T/ANICHZ6vD7WTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_spn(spn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b86e74e2-eaee-4393-aeeb-ae453653d230",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "test_data = np.array([1.0, 0.0, 1.0]).reshape(-1, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a6ef1b21-e0bd-4dbe-a026-fdafce049529",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.90730501]] [[0.14848]]\n"
     ]
    }
   ],
   "source": [
    "from spn.algorithms.Inference import log_likelihood\n",
    "\n",
    "ll = log_likelihood(spn, test_data)\n",
    "print(ll, np.exp(ll))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3a057fa3-b1a4-4c1a-af96-fe1e4217dcea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.68416146]] [[0.1856]]\n"
     ]
    }
   ],
   "source": [
    "llm = log_likelihood(spn_marg, test_data)\n",
    "print(llm, np.exp(llm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5805138c-e403-4dfc-b201-31fcb3c7ffef",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.68416146]] [[0.1856]]\n"
     ]
    }
   ],
   "source": [
    "test_data2 = np.array([np.nan, 0.0, 1.0]).reshape(-1, 3)\n",
    "llom =  log_likelihood(spn, test_data2)\n",
    "print(llom, np.exp(llom))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fadcdc48-fcee-408b-9538-74e73da3fab6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d5fabea3-4668-4359-ada5-2052e1869c42",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 1. 0.]\n",
      " [1. 1. 0.]\n",
      " [1. 1. 0.]\n",
      " [1. 1. 1.]\n",
      " [1. 1. 0.]]\n"
     ]
    }
   ],
   "source": [
    "from numpy.random.mtrand import RandomState\n",
    "from spn.algorithms.Sampling import sample_instances\n",
    "print(sample_instances(spn, np.array([np.nan, np.nan, np.nan] * 5).reshape(-1, 3), RandomState(123)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "4f990e8f-8c9e-4932-b590-8813ccd41127",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0.]\n",
      " [1. 0. 0.]\n",
      " [0. 0. 0.]\n",
      " [1. 0. 0.]\n",
      " [1. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "from numpy.random.mtrand import RandomState\n",
    "from spn.algorithms.Sampling import sample_instances\n",
    "print(sample_instances(spn, np.array([np.nan, 0, 0] * 5).reshape(-1, 3), RandomState(123)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f0853e35-5868-412a-8bef-d5f2a98da4d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "train_data = np.c_[np.r_[np.random.normal(5, 1, (500, 2)), np.random.normal(10, 1, (500, 2))],\n",
    "                   np.r_[np.zeros((500, 1)), np.ones((500, 1))]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "36c49359-5ab9-43c6-9a49-86697b975320",
   "metadata": {},
   "outputs": [],
   "source": [
    "from spn.algorithms.LearningWrappers import learn_parametric, learn_classifier\n",
    "from spn.structure.leaves.parametric.Parametric import Categorical, Gaussian\n",
    "from spn.structure.Base import Context\n",
    "spn_classification = learn_classifier(train_data,\n",
    "                       Context(parametric_types=[Gaussian, Gaussian, Categorical]).add_domains(train_data),\n",
    "                       learn_parametric, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "9c3791a3-97f3-4a5a-9a4d-74d6d17bd900",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_classification = np.array([3.0, 4.0, np.nan, 12.0, 18.0, np.nan]).reshape(-1, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "dbb925ac-474f-446b-aebc-68afefd78556",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.,  4., nan],\n",
       "       [12., 18., nan]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "560bb5b8-fc5f-4c4c-a25c-5837ab995a44",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.  4.  0.]\n",
      " [12. 18.  1.]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/zhangshen/software/anaconda3/envs/PY38/lib/python3.8/site-packages/spn/structure/leaves/parametric/Inference.py:88: RuntimeWarning: divide by zero encountered in log\n",
      "  probs[idx_in] = np.array(np.log(node.p))[cat_data[~out_domain_ids]]\n"
     ]
    }
   ],
   "source": [
    "from spn.algorithms.MPE import mpe\n",
    "print(mpe(spn_classification, test_classification))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "597ff357-7d8f-4d7d-8b7d-75278247432f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "eba0acc0-239b-49aa-8325-b917824eae43",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.          2.         -4.83457168  1.16542832]\n",
      " [ 1.          1.          3.44593668 11.44593668]\n",
      " [ 0.          2.          2.91977284  8.91977284]\n",
      " [ 0.          1.         12.49976944 15.49976944]\n",
      " [ 0.          1.          5.72734582  8.72734582]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "np.random.seed(123)\n",
    "\n",
    "a = np.random.randint(2, size=1000).reshape(-1, 1)\n",
    "b = np.random.randint(3, size=1000).reshape(-1, 1)\n",
    "c = np.r_[np.random.normal(10, 5, (300, 1)), np.random.normal(20, 10, (700, 1))]\n",
    "d = 5 * a + 3 * b + c\n",
    "train_data = np.c_[a, b, c, d]\n",
    "print(train_data[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "dd4d54ad-fb5c-49df-9bb1-eeebcd43881d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.    4.   -4.83   nan]\n",
      " [ 0.   12.   18.     nan]]\n"
     ]
    }
   ],
   "source": [
    "test_classification = np.array([3.0, 4.0, -4.83, np.nan, 0, 12.0, 18.0, np.nan]).reshape(-1, 4)\n",
    "print(test_classification)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "baeade35-d46b-4992-88d1-c6bf284e3490",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/zhangshen/software/anaconda3/envs/PY38/lib/python3.8/site-packages/spn/structure/Base.py:157: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
      "  self.domains = np.asanyarray(domain)\n",
      "/home/zhangshen/software/anaconda3/envs/PY38/lib/python3.8/site-packages/sklearn/cross_decomposition/_pls.py:97: ConvergenceWarning: Maximum number of iterations reached\n",
      "  warnings.warn('Maximum number of iterations reached',\n",
      "/home/zhangshen/software/anaconda3/envs/PY38/lib/python3.8/site-packages/sklearn/cross_decomposition/_pls.py:97: ConvergenceWarning: Maximum number of iterations reached\n",
      "  warnings.warn('Maximum number of iterations reached',\n"
     ]
    }
   ],
   "source": [
    "from spn.structure.Base import Context\n",
    "from spn.structure.StatisticalTypes import MetaType\n",
    "\n",
    "ds_context = Context(meta_types=[MetaType.DISCRETE, MetaType.DISCRETE, MetaType.REAL, MetaType.REAL])\n",
    "ds_context.add_domains(train_data)\n",
    "\n",
    "from spn.algorithms.LearningWrappers import learn_mspn\n",
    "\n",
    "mspn = learn_mspn(train_data, ds_context, min_instances_slice=20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "5f198d08-3aa9-40f7-bb40-5aa95fc7ca38",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.          4.         -4.83        0.76233681]\n",
      " [ 0.         12.         18.         24.35602634]]\n"
     ]
    }
   ],
   "source": [
    "# print(log_likelihood(mspn, train_data))\n",
    "print(mpe(mspn, test_classification))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "7bcb7a67-013d-458b-a377-ea3ba25214ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/zhangshen/software/anaconda3/envs/PY38/lib/python3.8/site-packages/spn/structure/Base.py:157: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n",
      "  self.domains = np.asanyarray(domain)\n",
      "/home/zhangshen/software/anaconda3/envs/PY38/lib/python3.8/site-packages/sklearn/cross_decomposition/_pls.py:97: ConvergenceWarning: Maximum number of iterations reached\n",
      "  warnings.warn('Maximum number of iterations reached',\n",
      "/home/zhangshen/software/anaconda3/envs/PY38/lib/python3.8/site-packages/sklearn/cross_decomposition/_pls.py:97: ConvergenceWarning: Maximum number of iterations reached\n",
      "  warnings.warn('Maximum number of iterations reached',\n"
     ]
    }
   ],
   "source": [
    "from spn.structure.Base import Context\n",
    "from spn.structure.leaves.parametric.Parametric import Categorical, Gaussian\n",
    "\n",
    "ds_context = Context(parametric_types=[Categorical, Categorical, Gaussian, Gaussian]).add_domains(train_data)\n",
    "\n",
    "from spn.algorithms.LearningWrappers import learn_parametric\n",
    "\n",
    "spn = learn_parametric(train_data, ds_context, min_instances_slice=20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ec703d7b-3d58-4876-a444-8a7d867a4bb7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.          4.         -4.83       -1.73256695]\n",
      " [ 0.         12.         18.         23.80580492]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/zhangshen/software/anaconda3/envs/PY38/lib/python3.8/site-packages/spn/structure/leaves/parametric/Inference.py:88: RuntimeWarning: divide by zero encountered in log\n",
      "  probs[idx_in] = np.array(np.log(node.p))[cat_data[~out_domain_ids]]\n"
     ]
    }
   ],
   "source": [
    "print(mpe(spn, test_classification))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "cb698269-e456-4c42-9ccc-35eb428c534a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-4.83914660187083\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "np.random.seed(123)\n",
    "\n",
    "train_data = np.random.binomial(1, [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.1], size=(100,10))\n",
    "\n",
    "from spn.structure.leaves.cltree.CLTree import create_cltree_leaf\n",
    "from spn.structure.Base import Context\n",
    "from spn.structure.leaves.parametric.Parametric import Bernoulli\n",
    "from spn.algorithms.LearningWrappers import learn_parametric\n",
    "from spn.algorithms.Inference import log_likelihood\n",
    "\n",
    "ds_context = Context(parametric_types=[Bernoulli,Bernoulli,Bernoulli,Bernoulli,\n",
    "                                       Bernoulli,Bernoulli,Bernoulli,Bernoulli,\n",
    "                                       Bernoulli,Bernoulli]).add_domains(train_data)\n",
    "\n",
    "spn = learn_parametric(train_data, \n",
    "                       ds_context, \n",
    "                       min_instances_slice=20, \n",
    "                       min_features_slice=1, \n",
    "                       multivariate_leaf=True, \n",
    "                       leaves=create_cltree_leaf)\n",
    "\n",
    "ll = log_likelihood(spn, train_data)\n",
    "print(np.mean(ll))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26911dc6-4d83-4f0e-8eb3-6411c42ba289",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "1bbc1df4-54ec-48fe-90e9-1110182c6f23",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Naive mle conditioning -4.326837245466242\n",
      "Random conditioning -4.309102071933407\n",
      "[[0. 0. 0. 0. 1. 1. 0. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 1. 1. 1. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 1. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 1. 1. 1. 0. 1. 0.]\n",
      " [0. 1. 0. 0. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 1. 1. 0. 1. 1. 0.]\n",
      " [0. 1. 0. 1. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 0. 1. 1. 1. 0.]\n",
      " [0. 1. 0. 1. 1. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 1. 1. 1. 1. 1.]\n",
      " [0. 0. 0. 0. 1. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 0. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 1. 1. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 0. 1. 0. 1. 0.]\n",
      " [0. 0. 1. 1. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 1. 1. 0. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 1. 0. 1. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 1. 1. 1. 0. 1. 0.]\n",
      " [0. 0. 1. 0. 0. 1. 0. 1. 1. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 1. 0. 1. 0.]\n",
      " [0. 0. 1. 0. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 0. 0. 1. 1.]\n",
      " [0. 1. 0. 1. 0. 0. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 1. 0. 1. 0.]\n",
      " [0. 0. 1. 0. 0. 1. 1. 0. 1. 1.]\n",
      " [0. 1. 0. 1. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 1. 0. 1. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 1. 0. 0. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 1. 1. 0. 1. 1. 0.]\n",
      " [0. 1. 0. 1. 1. 1. 1. 0. 1. 0.]\n",
      " [0. 0. 1. 0. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 1. 0. 1. 1. 0. 1. 1. 0.]\n",
      " [0. 0. 1. 1. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 1. 1. 0. 0. 0. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 0. 1. 1. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 0. 1. 1. 0.]\n",
      " [0. 1. 0. 1. 1. 1. 1. 1. 1. 1.]\n",
      " [0. 0. 0. 1. 1. 1. 0. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 1. 1. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 1. 1. 0. 0. 0. 0. 0. 1. 0.]\n",
      " [0. 0. 1. 0. 1. 1. 1. 1. 0. 0.]\n",
      " [0. 1. 1. 0. 0. 1. 1. 0. 1. 0.]\n",
      " [0. 0. 1. 0. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 1. 0. 0. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 1. 0. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 1. 0. 1. 1.]\n",
      " [0. 1. 1. 1. 1. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 1. 1. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 1. 1. 0. 0. 0.]\n",
      " [0. 0. 0. 1. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 1. 0. 1. 0.]\n",
      " [0. 1. 1. 1. 0. 1. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 1. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 1. 0. 0. 1. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 1. 0. 0. 1. 0. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 0. 1. 1. 0.]\n",
      " [0. 0. 1. 0. 1. 1. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 0. 1. 1. 1. 0.]\n",
      " [0. 1. 1. 0. 1. 0. 1. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 0. 0. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 1. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 1. 1. 1. 0. 1. 0.]\n",
      " [0. 0. 0. 1. 1. 1. 1. 1. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 1. 1. 1. 1. 0. 1. 0.]\n",
      " [0. 0. 0. 1. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 1. 1. 0. 0. 1. 1. 1. 1.]\n",
      " [0. 1. 0. 0. 1. 1. 1. 1. 0. 0.]\n",
      " [0. 0. 0. 0. 1. 0. 1. 1. 1. 0.]\n",
      " [0. 1. 1. 0. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 1. 1. 1. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 0. 0. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 1. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 1. 0. 0. 1. 0.]\n",
      " [0. 0. 1. 0. 1. 0. 0. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 1. 1. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 0. 1. 1. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 1. 1. 1. 0.]\n",
      " [0. 0. 1. 0. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 1. 1. 1. 1. 0.]\n",
      " [0. 0. 0. 0. 1. 1. 1. 1. 0. 0.]\n",
      " [0. 0. 0. 0. 1. 1. 1. 1. 1. 0.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "np.random.seed(123)\n",
    "\n",
    "\n",
    "from spn.structure.leaves.cltree.CLTree import create_cltree_leaf\n",
    "from spn.structure.Base import Context\n",
    "from spn.structure.leaves.parametric.Parametric import Bernoulli\n",
    "from spn.algorithms.LearningWrappers import learn_parametric, learn_cnet\n",
    "from spn.algorithms.Inference import log_likelihood\n",
    "\n",
    "train_data = np.random.binomial(1, [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.1], size=(100,10))\n",
    "\n",
    "ds_context = Context(parametric_types=[Bernoulli,Bernoulli,Bernoulli,Bernoulli,\n",
    "                                       Bernoulli,Bernoulli,Bernoulli,Bernoulli,\n",
    "                                       Bernoulli,Bernoulli]).add_domains(train_data)\n",
    "\n",
    "# learning a CNet with a naive mle conditioning\n",
    "cnet_naive_mle = learn_cnet(train_data, \n",
    "                            ds_context, \n",
    "                            cond=\"naive_mle\", \n",
    "                            min_instances_slice=20, \n",
    "                            min_features_slice=1)\n",
    "\n",
    "# learning a CNet with random conditioning\n",
    "cnet_random = learn_cnet(train_data, \n",
    "                         ds_context, \n",
    "                         cond=\"random\", \n",
    "                         min_instances_slice=20, \n",
    "                         min_features_slice=1)\n",
    "\n",
    "ll = log_likelihood(cnet_naive_mle, train_data)\n",
    "print(\"Naive mle conditioning\", np.mean(ll))\n",
    "\n",
    "ll = log_likelihood(cnet_random, train_data)\n",
    "print(\"Random conditioning\", np.mean(ll))\n",
    "\n",
    "# computing exact MPE\n",
    "from spn.algorithms.MPE import mpe\n",
    "train_data_mpe = train_data.astype(float)\n",
    "train_data_mpe[:,0] = np.nan\n",
    "print(mpe(cnet_random, train_data_mpe)) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b36ea3a-2ac2-4988-8797-ac0ea1d68d7e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
